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AN OPERATIONAL UNIFICATION OF FINITE DIFFERENCE METHODS 


FOR THE NUMERICAL INTEGRATION OF ORDINARY 


DIFFERENTIAL EQUATIONS 

By Harvard Lomax 
Ames Research Center 


SUMMARY 


One purpose of this report is to present a mathematical procedure which 
can be used to study and compare various numerical methods for integrating 
ordinary differential equations. This procedure is relatively simple, mathe- 
matically rigorous, and of such a nature that matters of interest in digital 
computations, such as machine memory and running time, can be weighed against 
the accuracy and stability provided by the method under consideration. 
Briefly, the procedure is as follows: 

(1) Find a single differential equation that is sufficiently represen- 
tative (this is fully defined in the report) of an arbitrary number 
of nonhomogeneous, linear, ordinary differential equations with 
constant coefficients . 

(2) Solve this differential equation exactly. 

(3) Choose any given numerical method, use it — in its entirety — to 
reduce the differential equation to difference equations, and, by 
means of operational techniques, solve the latter exactly. 

(4) Study and compare the results of (2) and (3) • 

Conceptually there is nothing new in this procedure, but the particular 
development presented in this report does not appear to have been carried out 
before. 

Another purpose is to use the procedure just described to analyze a 
variety of numerical methods, ranging from classical, predictor -corrector 
systems to Runge-Kutta techniques and including various combinations of the 
two. 


INTRODUCTION 


At present a large body of literature is devoted to the development and 
presentation of methods for integrating ordinary differential equations with 
given initial conditions. These methods are based on local polynomial 
approximation and are commonly divided into two classes, predictor -corrector 
methods and Runge-Kutta methods. The former are, as generally presented, not 



self -starting and use a fixed interval, or step, at which the function and 
its derivative are evaluated as the integration proceeds. The latter are 
self -starting and the interval of evaluation may vary from step to step. A 
current trend is to combine these two classes. The resulting methods are 
variously referred to as hybrid, generalized predictor -corrector, and 
combined. The latter designation is used herein. 

In this report a mathematical procedure, outlined in the summary, is 
presented which provides us with the capability of comparing these methods as 
they apply to simultaneous, linear, ordinary differential equations with 
constant coefficients. It is quite true that linear equations with constant 
coefficients are an extremely special set of all possible differential equa- 
tions, and, in fact, the numerical methods being discussed here are rarely 
used to solve them. However, such equations can be solved analytically both 
as differential equations, and as difference equations when transformed to 
the latter by a linear numerical scheme. The conclusion regarding the accu- 
racy and stability of a numerical method when studied in this way is, there- 
fore, precise. We need then only to defend the reasonable hypothesis that a 
numerical method which, on some given basis, is unquestionably inferior in 
solving linear cases, is, on the same basis, also inferior, in general , for 
use in solving nonlinear ones. 

When studied by the above procedure, all polynomial methods (known to 
the author) proposed for integrating ordinary differential equations fall 
into a smoothly connected system. By "smoothly connected," we mean, for 
example, that there is no sharp dividing line between predictor -corrector and 
Runge-Kutta methods. In fact, the standard, fourth-order, Runge-Kutta method 
is, in predictor -corrector terminology, a method composed of the successive 
application of an Euler predictor, an Euler corrector, a Nystrom predictor, 
and a Milne corrector. As such statements indicate, one of the principal dif- 
ficulties that can arise when different schools of thought are brought 
together Is the construction of a consistent and precise terminology. And 
the most troublesome problem in this area is to guard against conclusions 
based on implication. In particular, such a difficulty arises in the use of 
the term "step number" when combined methods are discussed. This is examined 
in the next paragraph. 

All numerical methods of the type being considered are cyclic in appli- 
cation; that is, for a fixed reference value of the independent variable, 
a pattern of calculations is performed (solving equations, evaluating deriv- 
atives, estimating errors, etc.). At the end of these calculations the 
value of the function has been determined at a point advanced by some inter- 
val. The independent variable is re -referenced ahead by another interval and 
the identical pattern of calculations is repeated. These cycles are con- 
tinued indefinitely . The interval Involved is referred to as the step size. 

The number of locations, spaced by a fixed value of this interval , at which 
the function and/or its derivative are retained for use in the next cycle or 
pattern of computations, corresponds to the step number of the given method. 

The definitions of step size and step number when made in this way hold for 
these terms as they are generally used in the literature for both combined 
and uncombined methods. In predictor -corrector schemes this step number is a 
fundamental parameter that can be, and is, used to connect the stability and 
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accuracy of a given method. In fact, in a well known theorem, Dahlquist 
states that no stable , predictor - corrector method with a step number, k, can 
give a local polynomial approximation of order k+ 2 if k is odd, or of 
order k + 3 iY k is even. However, in Runge-Kutta methods, or methods that 
combine the Runge-Kutta and predictor -corrector concepts, this step number is 
not connected in any way with stability . Thus, stable combined methods having 
any value for the step number (including one) can be constructed that will 
fit local polynomials of any order. This implies that combined methods are 
greatly superior to uncombined ones. But, in fact, for fixed values of 
machine memory and running time, the maximum order 1 of a local polynomial fit 
appears to be the same for stable methods combined or uncombined. 

At the beginning of the report certain basic terms are defined so as to 
make the subsequent discussion more precise. Then the approach to be used in 
the analysis is described and it is shown that a single representative dif- 
ferential equation can be used to study the accuracy and stability of 
difference-differential approximations as they apply to the analysis of simul- 
taneous differential equations. An attempt is made to classify various class- 
ical and modern numerical methods according to three categories: 

1) The number of iterations per cycle of computation 

2) Whether they are complete or incomplete 

3) Whether they are combined or uncombined. 

Some general procedures falling into certain combinations of these categories 
are analyzed in detail. Finally, the operational form of a difference- 
differential equation is defined and its implications with regard to the 
study of numerical methods is discussed. 


SYMBOLS 


A constant in representative equation (See eq. (37)-) 

DE(E) see equation (52) 

erp error of a numerical method in terms of local polynomial approxi- 

mation 

er^ error of a numerical method in calculating the particular solution 

of the representative equation (37) (See eqs. ( 57 ) and ( 63 ).) 

er^ error of a numerical method in calculating the complementary solu- 

tion of the representative equation (37) (See eqs. ( 68 ) and 

(71)-) 


1 The magnitude of the leading error term found by means of a Taylor 
series expansion is lowest for the combined methods, however. 
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difference operator (See eq. (9)0 

computational step size (See eq. ( 123 ) 0 

representational step size (See definition (2)0 

index used in expressing difference -differential equations 
(See eq. (3) • ) 

k + 1 - j 

step number in a predictor or corrector 
coefficients in the operational form (See, e.g., eqs . 

(51)0 

reference step location 
see equation (52) 

coefficients in the operational form (See, e.g-, 

eqs. (51)0 

dependent variables 

du dvr 
dx 9 dx 

independent variable 

coefficients of dependent variable in difference- 
differential equations 

coefficients of derivative of dependent variable in 
differential equations 

representative eigenvalue, that is, coefficient of u in 
representative differential equation (37) 

spurious roots of difference equation 

principal root of difference equation 

induced stability boundary (See eq. (73)0 

representative maximum frequency (See eq. (37)0 

eigenvalues of simultaneous ordinary differential 
equations (ll) 
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DEFINITION OF TERMS 


Some of the following expressions are in common usage but vary slightly 
in meaning with different authors. The definitions given below are intended 
for this report to simplify and make more precise the subsequent discussion. 

Difference-differential equations : 

Let the dependent variable u be a 
function of the independent variable 
x. Let u T represent the derivative 
of u with respect to x and desig- 
nate x n by nh and u(x n ) by u n 
where n is an integer and h is a 
constant. Then equations which 
relate u n+ ^ +1 _ j , un+k+i - j and 

x n+k+i - j where j = 1, 2, . . 
k + 1 are called difference- 
differential equations with step 
number k . 

Predictor : Any difference-differential equation relating u n +k to 

values of u and u T at previous steps. Thus, for a k-step predictor 

u n+k ” ^( u n+k+i-jj u n+k+i-j^ x n+k+i-j); 3 ” 2 , • • •> k + 1 

A predictor is an explicit formula that extrapolates given data. 

Corrector : Any difference-differential equation relating u^^ to the 

values of u and u T at n + k as well as to those at previous steps. 
Thus, for a k-step corrector 



*n 

Sketch (a) 


u n+k “ f ( u n+k+i-j j u n+k+i-j> x n+k+i-j) } j - 1, 2, . . . , k + 1 

In this form the corrector is an implicit formula. In practice the 
values of u and u r in the arguments • of f are generally those 
determined by predictors or previous correctors. 


Iteration : In the numerical solution of ordinary differential equations 

the repeated calculation of the right-hand side of equations having the 
form 

u ! = F(x, u) (la) 


or for multiple equations 

t 

Ui 

T 

u 2 


n(x, uu u 2 , 
F 2 (x, u 1# u 2 , 



(lb) 
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is necessitated. In this report we refer to every such evaluation (i.e., 
explicit calculation of the derivatives using the differential equations) 
as an iteration. By this definition, methods composed only of predictors 
require one iteration per step. Methods using one predictor followed by 
one corrector require two iterations per step, etc. 

Reference step : We will inspect a wide variety of methods in which the 

words "step size," by common usage, have different implications. In 
order to have a parameter by means of which all methods can be compared 
on a common basis, the term "reference step" is introduced and designated 
by the symbol H. 


H = the increment in x that a solution /px 

Is advanced by two iterations ' 

In many applications the numerical calculations necessary to evaluate 
the derivatives, Fj(x,u) in equations (l), are extremely complicated and 
time consuming. In such applications, if errors are referenced to H, 
the accuracy of various numerical methods can be compared with the assur- 
ance that the total machine running time will be very nearly the same. 
Since most methods in practical use employ a predictor followed by just 
one corrector, two iterations were chosen for a base (rather than one) 
so that H would coincide with the most commonly used error reference. 
Both 2 h, the computational step size, and H, the reference step size, 
are used in error terms in the following analysis. 


Cycle of computation : All the calculations and logic required to advance 

the data while n refers to the same location. A cycle is completed 
when all the dependent variables and their derivatives at n + k have 
been calculated as accurately as the chosen method permits and 
preparation for stepping ahead commences. 


Family : Any combination of values of u, u T and other families at 

n + k, n + k - 1, - . . , n that is formulated and used in a cycle of 
computation. A family may or may not be saved for future cycles of com- 
putation. In this report a family is usually designated by a super- 
script, and a predictor always generates the first family. A derivative 
belongs to that family of u used in its calculation; that is, 



F 



) 


Final family : The new values of u and u* last evaluated in a cycle 

of computation. The superscript is always omitted from the final family 
of u (its distinguishing feature) and sometimes from the final family 
of u r (see the definitions below of complete and incomplete methods). 


^When corrector methods are analyzed as such, without regard to how their 
equality is brought about, the reference step H is undefined. 
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Memory in a k-step method : All those values of u, U* and other 
families (if there are any) that are used but do not change during one 
cycle of computation. 

Incomplete methods : Methods for which the dependent variable and its 

derivative are members of the same final family. That is, after the 
dependent variable is evaluated for the last time at a given point, it 
is used to calculate the derivative at the same point. Most "conven- 
tional” methods (Hammings 1 s, Milne* s, etc.) are incomplete. In this 
case the superscript is omitted from the final family representing the 
der ivative. 

Complete methods : Methods in which the derivative of at least one final 

family is never evaluated. They are referred to as complete because 
they most completely fill the matrix which determines the operational 
form. 

Combined methods : Methods that combine the concepts usually separately 

designated as predictor -corrector and Runge-Kutta. A combined method 
can be thought of either as a predictor -corrector method without equal 
spacing, or a Runge-Kutta method with memory. Combined methods can be 
either complete or incomplete. 

Fundamental family : One that is computed using a memory composed only 

of final families. 

Embedded polynomial : The highest order polynomial which is an exact 

solution to a given set of difference-differential equations. 


FtMDAMENTALS 


Difference-Differential Equations 
Two of the simplest difference -differential equations are 

u n +i - ^ - hu^ = 0 

and 

u n+i " % " \ h ( u A+i + u n) ^ 0 

and are referred to as the Euler and modified Euler equations, respectively. 
These and all such formulas presented in books on numerical analysis are 
special forms of the general, linearized, k-step, difference -differential 
equation with constant coefficients which can be written 

u n+k - Pi u n+k " kPi u n+k - • • • Pj u n+k-j ” ^Pj u n+k-j “ • * * 


- Pk+i u n - h 3k+i u A = 0 (3) 
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Nearly always, equation (3) represents formulas based on polynomial 
approximation. This simply means that if each u and u r is expanded in a 
Taylor series about x n , the coefficients of the powers of h in equation (3) 
will vanish up through some integer L. The number L is then the order of 
the polynomial approximating the function in the interval x n < x < x n+ k 
(the embedded polynomial) and the product of h^ 1 , and its coefficient is the 
first term of the truncation error. For example, since 

u n+1 = u(nh + h) = u(x n + h) = u n + hu^ + ^ h 2 u£ + g h 3 u£’ + . . . 

hu^ +1 = hu'(x n + h) = hu^ + h 2 u^ + ^ h 3 u^' + . . . 
we can construct for the modified Euler method the simple table: 


From 

U n 

hu£ 

h 2 u£ 

h 3 u^’ 


u n+i 

1 

1 

1/2 

1/6 


-u n 

-1 

0 

0 

0 

• • • 

-l/2 hu^+i 

0 

-1/2 

-1/2 

-i/4 

• • • 

-1/2 h< 

0 

-1/2 

0 

0 

• • • 

Sums to 

0 

0 

0 

-1/12 

| • • • 


Clearly, the order of the polynomial embedded in the modified Euler method is 
2 (even though only one step is used) and the truncation error is predomi- 
nantly -u ,M h 3 /l2. A similar tabulation for equation (3) is shown below. 


From 

U n 

huA 

h 2 4 

h 3 u^ 

* h 4 ufi' 

• 1 

u n+k 

1 

k 

1/2 k 2 

1/6 k 3 

1/24 k4 


-Pi u n+k 

"Pi 

-Pik 

-1/2 3 x k 2 

-1/6 Pik 3 

-1/24 Pik 4 



0 

-Pi 

-Pik 

-1/2 Pik 2 

-1/6 pik 3 


“02 u n+k-i 

-02 

-Pa(k-l) 

-1/2 P 2 (k- l) 2 

-1/6 p 2 (k- l) 3 

-1/24 p 2 (k- l) 4 


-hP^U-n+k-i 

0 

-Pa 

-P a(k - 1) 

-1/2 p 2 (k- l) 2 

-1/6 P4(k-1) 3 


-Pk+i u n 

-Pk+i 

0 

0 

0 

0 


-hPk+i u n 

0 

-3k+i 

0 

0 

0 


Sums to 

0 

0 

0 

0 

0 

. . . 


Equating the sum of the first L columns to zero gives the conditions on the 
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P.? and p! required if equation (3) is to represent a polynomial of order L 
through J k + 1 points. One can show that the product of the sum of the Zth 
column and its heading is 


er p ( l ) -- 



k+x 






1 - jf'Uj + Pj(k + 1 - j) Z 



(4) 


Therefore, er^(L + l) is the first term in the truncation error of a Taylor 
series expansion for any function u(x) represented by the differ ence- 
differential equation (3) • The total truncation error is given by 

oo 

^ er p (z) 

Z=L+1 

Until recent years 3 equation (3) was used as the sole basis for deter- 
mining the accuracy and stability of a numerical method. As is now well rec- 
ognized, this is rarely a correct procedure. Let us suppose, for example, we 
are using equation (3) to find the value of u at x n+ ^. Then the only time 
it describes the total numerical method is when both and P x are identi- 
cally zero. This is the case when a predicted value is calculated but no 
correction is made. Equation (3) also represents the numerical result of an 
implicit method where the terms multiplying Pi and p{ might have been cal- 
culated by some iterative procedure. This is the assumption under which it 
is usually applied. Almost all practical methods use at least one corrector 
and when such is the case the accuracy and stability of the actual results 
are affected by the mutual interaction of the predictor and all of the sub- 
sequent correctors. This will be fully developed in the following sections. 

At this point we wish merely to define a notation for a true predictor - 
corrector process. Consider, as above, that the values of u n +k and u^^ 
are unknown but all values with prior subscripts are known, being either 
given or obtained from previous calculations. Then for the predicted value 
at n + k we can write 

( a j u n+k+i -j + a j hu n+k+i-j) (5) 


k+i 

( a j u n+J + a jlxun+j) (6a) 


(x) 

u n+k = 


k+x 




or , to shorten the notation. 


(i) 

u n+k ~ 


See the next section for some 


historical discussion. 
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where J = k + 1 - j. This forms the first family 4 at the location n + k 
which is designated by the superscript (l) . If this is followed by a cor- 
rected value at n + k, we can write, where again J = k + 1 - j. 


(2) 

%+k 


k+i 

“ Pi hu nik + Y ^j Un+J + 


(6b) 


J- 2 


which defines a second family at n + k. (The possibility of including a 3i 

term is discussed under equation (44).) Another corrector could be added with 

coefficients 7 . 8 , 7 ! forming a third family, etc. If, however, we consider 
J 3 

the cycle of computation complete after the evaluation of equation (6b), then 
the second family is the final family for the dependent variable u. Next, 

f 2 ) 

a decision must be made as to whether or not u^ + k shall be used to evaluate 

another estimate of the derivative at n + k. (in this regard, see 
Definition - Family.) A derivative has already been calculated at n + k, 

(1) 1 

namely u^+k , and it can be used to advance the solution* If this path is 
chosen, we (in this report) refer to the method as a complete method. The 
function and its derivative are members of different families (initially 
generating the various families presents the same difficulties that arise in 
starting multistep methods) and equations (6) can be written 

k+i 


u. 


h) 

ri+k 


I (“j 


1 j u n+ J 




<3 ="2 
k+i 


u n+k 


p i4+k 


Y ( P J Un+J + P >ntr') 


J- 2 


(7a) 


where the superscript (2) has been omitted from the final family for u. 
Choosing the other path provides (what is referred to herein as) an incomplete 

(2) (2) ! 

method. In this case Un+fc is used to find and the latter is placed 

in memory for use by subsequent predictors and correctors. Now, the function 
and its derivative are members of the same family and the superscript (2) is 
omitted from both; thus, 


%+k - 


d t ( 1 ) 1 
u n+k = Pi h Vk 


k+i 

Y ( a j u n+J + a j hu n+j) 


k+i 


(^"n+j + ^ hu n+j) 


J=2 


(7b) 


Incomplete methods are the most common, but not necessarily the best. 
4 See Definition of Terms. 


10 



Table I lists the coefficients of a few of the commonly used difference- 
differential equations, together with the leading terms of their truncation 
errors as defined by equation (k) . Identifying names axe included, although 
these are not unique. 


The Method Used to Measure Accuracy and Stability 

Development . - In the application of equations (7) to equations (l) , the 
matter of overall accuracy in the resulting numerical scheme depends not only 
on the truncation error but also on the stability of the numerical process as 
it is continued along a number of steps. Thus, as is well known, the modified 
Euler method, row 1 in table l(b) , is stable; but the Nystrom equation, row 2 
in table l(a), when used by itself, is unstable. The usefulness of any set 
of difference-differential equations depends upon a balance between accuracy 
and stability. Dahlquist (ref. l) found the maximum order of a polynomial 
that could be embedded in equation ( 3 ) for a given k under the condition 
that the resulting predictor -corrector method would be stable as h -> 0. He 
concluded, for example, that a three-step method could never support a poly- 
nomial of order h 5 or higher and still be stable. Hamming (ref. 2), at 
about the same time, developed a stable three -step corrector formula having a 
truncation error led by a term of order h 5 , the minimum possible according 
to the proof of Dahlquist. Hamming’s stable corrector formula and the most 
accurate, but highly unstable, three-step corrector formula are shown in rows 
6 and 9 in table l(b) . 

In a very interesting paper, Chase (ref. 3) put a new light on the 
developments mentioned above and brought out two important points. First 
(with a notable exception in Hamming’s article) nearly all theorems and proofs 
regarding stability published prior to Chase’s paper are based on the limiting 
case when h -> 0. Second, nearly all analyses of corrector formulas, includ- 
ing Hamming’s, assume that all effects of the predictor have vanished, or, in 
other words, that the corrector equation is brought into complete balance. 
Chase showed that when the above conditions are not met (which is certainly 
the practical case, since step sizes cannot be zero and very often the cor- 
rector is used only once) , the conclusions regarding stability of the various 
methods undergo startling changes. For example, Hamming’s corrector formula - 
fully satisfied - is stable for values of -|7vh| even less than -2. But 
the result of predicting with row 6 of table 1 (a) and correcting with row 6 
of table l(b) (Hamming’s method without the modifiers) is unstable for 
— | Ah | < -0.5- Even more dramatically, Chase showed that predicting with row 
6 in table l(a) and correcting with row 5 of table l(b) (Milne’s method with- 
out modifiers) is stable for -0-3 > -|Xh| > - 0 . 8 . In other words, a predic- 
tor and a corrector individually unstable for all — | Ah | <0 combine to form 
a stable method for a certain applicable range of - | Ah | ! 

A final lesson to be learned from Chase’s paper concerns the use of 
’’modifiers.” These are weighted combinations of the predicted and succesive 
corrected values. In the two cases mentioned above, Hamming’s and Milne’s, 
the use of modifiers in one case increased and in the other case decreased 
the range of stability. This suggested a study which would determine optimum 
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values for the modifiers. In the terminology of the present report, a modi- 
fier is a family that is not fundamental. By means of operational techniques, 
we shall see that a method with a given step number and modifiers can be iden- 
tified with a method having a higher step number but composed only of 
fundamental families . 

Description . - Consideration of the matters discussed above suggested an 
approach differing from the one presented in references 4 and 5 . The 
approach in this report does not optimize the coefficients p* and P 1 in 

equation (3), but rather the coefficients in an "operational form" as defined 
in the last section of this report. The coefficients defined are the funda- 
mental parameters governing the accuracy and stability of linear numerical 
quadrature formulas and are functions of all the aj, aj, P j, and pj and their 
complete interaction. Furthermore, stability will not be studied in the limit 
as h -*• 0, but rather over a finite range of h, a range which is to be made 
as large .as possible for a given accuracy, and includes the entire complex 
plane . 

Sketch (b) graphically presents the details of the proposed approach. 

A representative differential equation (or set of differential equations), 



which is discussed fully in the next section, is chosen and solved exactly. 
Then a group of linear difference -differential equations with unspecified 
constant coefficients is introduced and combined with the differential equa- 
tion to form a set of linear difference equations. Operational techniques 
are used to solve these difference equations exactly. The solution to the 
difference equations is then compared with the solution to the differential 
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equations. Eventually, the coefficients y , etc., are chosen so that 
the two exact solutions match as closely as possible under the condition that 
the resulting process remain stable over a given range of step size. 

Operational Solution of Difference Equations 

In the presentation of the analysis in the following sections the reader 
is assumed to have some knowledge of the theory of ordinary difference equa- 
tions. This theory is well developed but its publication is not nearly so 
widespread as that on the theory of differential equations. A brief review of 
a portion of an operational approach to difference equations is given below. 
For complete treatments see Boole (ref. 6) or Milne-Thomson (ref. 7 ) • 

Classically, 


u n +k + C!U n+ k_i + . . . + C k u n = F(n) (8) 

is defined as an ordinary difference equation of order k. Unfortunately, 
the word order in most modern books and articles on numerical methods is used 
to designate the highest integer exponent in the polynomial embedded in equa- 
tion (8) . In this report we also refer to the order of a method in the latter 
fashion and refer to k in equation (8) as the step number. If the coeffi- 
cients Ci, C 2 ; . . . , C^ in equation (8) are independent of n and u, the 
equation is an ordinary linear difference equation of step number k with 
constant coefficients. The solution to the equation when F(n) = 0 is the 
complementary solution which is added to the particular solution involving 
F(n) when F(n) ^ 0. 

Solutions to equation (8), when the coefficients are constant, are 
obtained by operational methods similar to the Laplace or Fourier transforms 
for differential equations. If E is defined as the operator 

Eu n = e h ^ d / dx) u n = u n+1 (9) 

equation (8) can be written 

[E k + CnE 11 " 1 + . . . + C k }u n = F(n) 

The complementary solution is determined by finding the roots to the charac- 
teristic equation 

E k + CxE k_1 + . . . + C k = {E - Ax}{E - A 2 } . • . {E - A k } = 0 
and this solution is 

— _n -=• ..n 

Un — C x Ai + C 2 A 2 + . • • + C^-A^ 

where the Cj are constants determined by the initial conditions. If a root 
is repeated m times, its coefficient is a polynomial in n of order m - 1. 
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Thus , for 


[E - A;L} m+1 {E - = 0 


the solution would be 

u n -(Cio + Cxin + • • • + + C2A2 

If we consider forms of equation ( 8 ) for which the Cj are real, any complex 
root of the characteristic equation must have a conjugate. Such cases can be 
treated as in the following example. Let 

[E - a - ip}{E - a + ip]{E - A 3 } =0 

be the characteristic equation. Then setting 

r 2 « a 2 + (3 s 


0 ~ tan(p/a) 


we can write the equation in the form 

u n = G x r n cos nO + C 2 r n sin nS + C3A3 


The particular solution to equation ( 8 ) is easily expressed when F(n) is 
a polynomial in n or any sum of terms having the form y n where y is any 
complex constant. The latter solution is given by Boole ! s first rule and is 
of particular value to us . Boole showed 

{G(E)}~ n = 7^(7) ( 10 ) 

where the notation reads G(E) operating on 7 11 equals the product of ~ n 
and G(~) • Thus, the particular solution of 

{E k + CiE k_1 + . . . + C k }u n = Ae^ hn 


is 


( u n)p 



+ CiE k-1 


1 

+ 



Ae 


qhn 


or 


/ x A M'hn /{ qhk 
(u n )r) = Ae /(e 




+ c k ) 


It is valid for complex [x and its evaluation does not require that the 
roots to the characteristic equation be known. 

The solution of simultaneous linear difference equations with constant 
coefficients offers no particular difficulty. For example, the two equations 


111 - 



u n+i + 2v n = 2 


-n 


u n+i - \ u n + v n+i ~ 0 

in operational form become 

Eu^ + 2v n = 2 
{ E “ "n + Ev n = 0 


or in matrix notation 



The characteristic equation for the set is 


E 


E - | E 


= E 2 - 2E + 1 = (E - l) 2 = 0 


and the general solution for u n is 


u n “ (^o + Ci n ) ( + l) + 


E 


or 


— Cq + Gin + 2 


(E - 1)‘ 

(i-n) 




THE REPRESENTATIVE DIFIERENTIAL EQUATIONS 


Development 

Fundamentally, the representative differential equations, for which the 
conclusions throughout this report exactly apply, are the set of simultaneous, 
linear, first-order, differential equations with constant coefficients 
( w T = dw/dx) 
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“ a n w i + ai 2 w 2 

T 

W ^ — ^ 21^1 + a 22^2 
w 3 - a 3l v l + a 32 w 2 


+ fl(x) 

+ f 2 (x) 

+ fs(x) > 


(11) 


or any group of higher order differential equations which reduce to such 
a set. Although the subsequent analysis rigorously applies to equations (ll), 
it is unnecessary, for the purpose of studying the accuracy and stability of 
a numerical method, to consider them in such a general form. In fact, we will 
see in the next few sections that the stability and accuracy that result from 
integrating equations (ll) by any of the numerical methods considered herein 
are completely independent of the elements a-M except as these elements 
determine the roots of the characteristic equation (the eigenvalues) . In 
other words, if equations (ll) are integrated by some polynomial numerical 
method for any number of steps and then uncoupled (put in the form of eigen- 
vectors), or if equations (ll) are first uncoupled and then integrated using 
the same method and step location, the results will (except for roundoff 
error) be identical regardless of whether or not they are correct or the 
numerical method is stable. 

Exact solution using the Laplace transform. - To begin with, let us con- 
sider the exact solution to equations (ll) as it is derived by means of the 
Laplace transform. If w(s) is the Laplace transform of w(x) 


00 


w(s) = J e _sx w(x)dx 

o 


(12) 

and 



sw(s) - w(0) + f° e’ SX w ! (x)dx 
^0 


( 13 ) 

Multiplying both sides of equations (ll) by e“ sx and integrating with 
respect to x from 0 to ^ gives 


(s-n - s)w^ + ai 2 W2 + . • . = -W]_(0) 

a21 w l ( a 22 *” S ) W2 = -W2(0) 

3 * 31^1 & 23^2 ~ - W3 ( 0 ) 

- L(s) 

- L(s) 

- f* 3 (s) 
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a set of simultaneous, linear equations for Wj., w 2 , . . . w m , .... The 
determinant of these equations is a polynomial in s; and equating this poly- 
nomial to zero results in the characteristic equation for the differential 
equations (ll) . Let the roots to this characteristic equation be 
cfi, a 2 , • • • .... (For simplicity, the argument proceeds as if none 

of the roots is multiple, but this restriction is not at all necessary.) From 
the well -developed theory of the Laplace transform, the complementary solu- 
tion to the equations (ll) can at once be written 


n - c ii eCTlX + c al e ° sX + ■ 


+ ♦ 


or since x =• nh 


“ CizCe 01 ^) 11 + C^z ( e astl ) + . • . + C m 2(e 0ml1 ) + . . . (14) 

where C m i are constants depending on the initial conditions and the nature 
of the functions fi(x) . 

As is well known, the above developments can be viewed in a slightly 
different light. Write equations (ll) in the matrix form 

w T ~ [A] w + f (15) 

where -» defines a column vector and the [A] matrix is defined by 


9*11 

a i2 

ai3 

• • m ' 

a 21 

a 22 

a 23 

. . . 

a 3l 

a 32 

a 33 

... 


[A] = 


(16) 


The complementary solution (l4) immediately follows where the cr m are the 
eigenvalues of [A]. 

Numerical solution by a predictor or an implicit corrector . - Now let us 
solve equations~(ll) using a single difference-differential equation. In 
practice this would be a predictor or corrector with the implicit relation- 
ship somehow brought to equality. In the next section the generalization of 
the following to actual predictor -corrector methods will be discussed. The 
analysis is presented in this order because of the simplicity of the develop- 
ment in this section relative to that in the next. 

Consider the difference -differential equation (recall that J= k+ 1 - j) 

k-fi k+i 

u n+k =* ^ P.j u n+J + h ^ Pj u n+J ( 1 7) 
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which is a predictor if Pi - Pi - 0. Introduce the operator E (see eq. (9)) 
and rearrange. There results 


or 



(18) 


which defines the operator S in terms of the operator E. Thus, if the 
difference-differential equation is applied to the differential equations (ll), 
there results at the nth step the set of linear difference equations 

(an - S) wm + a i2 w 2n + . . - = -fi(nh) 
a 2i w in ( a 22 " S) w 2n ~ — f 2 (nh) 

a 31 wm + a 32 w 2 n - -f 3 (nh) 


Clearly, the roots to the characteristic equation in S are once again the 
eigenvalues of the matrix [A]. In other words 

(s - Ci ) (S - cr^) • • • (8 — o'm) • • • = 0 

This leads to the rather remarkable result that the numerical method isolates 
the individual roots of the exact solution and operates on each of them 
individually as_ if the others were not even present ! 

Recall the definition of S and construct the "subcharacteristic" 
equation for E in terms of a m . Thus 

k+i 

E k - ^ (3j + a m hpJ)E J - 0 

j - 1 

which has the "subroot" structure 

(E - Aim)(E - Apm) . • ■ - 0 
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One of these roots will approximate the Taylor series expansion of the term 

with a truncation error appropriate to the degree of the polynomial 
embedded in the difference -differential equation. This root is commonly 
referred to as the principal root, designated Amu and can always be 
expressed in the form 

Aim “ 1 + tfrnh + \ °mk 2 + \ ^m^ 3 + * • • - e m ± (19) 

o' h 

where the ± after the e° m indicates the existence of a truncation error. 
The remaining roots A 2 nu A 3m , • • • az’e spurious. They are introduced by the 
numerical method and depend on the choice of the difference -differential 
equation, both as to number and magnitude. 


In summary, the exact solution to the differential equations is 

"Zn = C lZ (e° lh ) n + • . • + C ral (e ffmh ) n + . . • 
and the exact solution to the difference equations is 

W ln - Cxi,t(e" lh ±) n + • • • + ♦ . . . 

+ C 2 i,z(A2i) n + . . . + C 2m ^2(A2m) n + • • • 

+ ^31, Z (A3i) n + • • . + C 3m ^(A 3 tn) n + • • • 


( 20 ) 


( 21 ) 


The results include complex roots and may be extended to multiple roots. 

Numerical solution using both a predictor and corrector .- The conclusions 
drawn in the previous section were shown to be true when a single difference- 
differential equation is used to integrate a set of simultaneous differential 
equations. We now show the same conclusions hold for a predictor -one - 
corrector method. 


To begin with, apply the incomplete predictor -corrector combination 


(i) _ 

u n+k “ 


k+i 


^ ( a j u n+J + a j hu n+j) 


j= 2 

T k+1 

u n+k = Pi hu n+k + ^ (PjU n+ j + pjhu^+j) 

J= 2 


( 22 ) 
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to the single equation 

u T - cr u 

There results, in operator notation. 



Define a new set of operators 


s - 



t 


piE k 


k+i 





(23) 


(24) 


and the matrix equation 


— 


(1) 

s a - y 


u n 

at cr - z 


1 


(25) 


follows. The characteristic equation is simply 


s(cr - z) - crt (cr - y) = 0 (26) 

Substituting for s, t, y, and z one can determine the root structure for 
E in terms of the eigenvalue of the single equation (23) • This root struc- 
ture determines the accuracy and stability of the method defined by 
equations (22) as it applies to a single differential equation. 
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For the two equations 


u i — 3 -ir u x + ^± 2 ^ 2 . 


Up — a*2x u x + a ^p Up 


(2 7) 


the matrix expands to 


s 

0 

a n - x 

a 12 


-con 

u-m 

0 

s 

a 21 

app — x 


(0 

u 2 n 

axit 


a n - y 

a 12 


Tixn 

a 2 it 

9 >Ppt 

a 21 

a 22 - y 


u 2 n 


= 0 


( 28 ) 


If the eigenvalues of [A], where 


[A] 


■[ 


a ll a 12 

a 21 a 22J 


are and cr £ , one can use the identities 

s ±8 2. — a ll a 22 “ a 21 a 12 


CTi + ffg — an + 8*12 

2 2 2 -, 2 
Cx + O' 2 - a XX + 2 axx a 21 + a 22 


and show that the determinant of the matrix in equation ( 28 ) reduces to the 
product 


2 s °j “ 7 

n 

j=l ffjt CT-j - z 


(29) 


having the characteristic equation 


[s(crx - z) - tffxCffx - y)][s(cr 2 - z) - ta 2 (cr 2 - y) ] =0 (30) 


This shows that when the method defined by equations (22) is applied to two 
coupled differential equations > the stability and accuracy of the result 
depend entirely on the eigenvalues of [A] and are not related (in any other 
way) to the magnitude of the individual elements . 

The generalization of the . above to larger groups of simultaneous equa- 
tions depends upon a proof of the conjecture 


21 


s 

0 

0 


a-n-x 

a 12 

a 13 



0 

s 

0 


a 21 

a 22 “ x 

a 23 

• . . 


0 

0 

s 


a 31 

a 32 

a 33 -X 



8-1 it 

a i2 t 

ai3t 


an-y 

a 12 

a 13 


“ n[s(dj-s)-tdj(d ;5 -y) ] 

a 2 it 

a 2 2t 

a 23 t 


a 21 

8-22-7 

a 23 

* * ‘ 


as it 

a3 2 t 

a33t 


8-31 

a 32 

a-33 -y 




' • ... , ( 31) 

where the Oj are the eigenvalues of the A matrix (see eq. (l6)). Although 
the author knows 5 of no formal proof of (31) , its equality can easily be 
demonstrated, within the bounds of numerical error, by means of a digital 
computer. Simply evaluate and compare the left- and right-hand sides of (31) 
(subroutines that calculate the eigenvalues of real matrices are standard 
equipment) using random numbers for the terms t,s,y,z and the elements ap-j. 
This "numerical proof" has been carried out for a variety of determinants 
with order 2 through 10. 


The preceding analysis leads to a final hypothesis. First consider the 
following : Let 



s Knew at this writing (see end of this section). 


O O 


( 32 ) 


where by and Ay are constants. Then if 

w’ — [Al* + [B]f 

where 

[A] = [B3[C][Br (33) 

a is independent of [B] where 

u = [B] 1 w ( 3 I 4 -) 

The proof follows immediately by multiplying equation (32) by [B]” 3 ". How 
reduce equation ( 32 ) to a set of difference equations by applying any com- 
bination of linear, difference-differential equations representing complete 
or incomplete, predictor, multiple-corrector methods. Solve these equations 
for U and apply equation (3^-) • Regardless of the choice of [B], we hypoth- 
esize that the numerical value of u will be identical at every step except 
for roundoff error. Again a ’’numerical proof” of this conjecture was obtained 
by solving five linear, simultaneous equations with real and complex eigen- 
values for several choices of [G] and f and a variety of [B]. Fith double- 
precision arithmetic the differences in u caused by the various choices of 
[B] could safely be attributed to roundoff. 


During the preparation of, this report a simple analytical proof of the 
above was presented to the author. This proof, which changes the hypothesis 
into a theorem, was prepared by Dr. Filliam A. Mersman 6 and proceeds as 
follows . 

Combine equation (15) with the predictor -corrector sequence defined by 
equations (22) . Using [i] to designate the unit matrix, we have 


Hn+tc 


+ hy [A])w n+J + haj[I]f u+ j| 


,( 1 ) 


3 

k+i 


k+x 


^n+k = Pxh[A]^ii + + hPj[A])^ n+J Uh ^ Pjf a+ J 

j=2l J j'i 

_Jx) 

Eliminating w n+ j and introducing an operational notation, we derive the 
equation 


k+i 


'^[I] - ^ (PjCU + h-( Plocj + Pj)[A] + h 2 p{aj[A] 2 )E J J^w, 

J = 2 

r k+x 

■ h u 


n 


s+x k+x A 

^ P][I]# + Pxh[A] ^ 

1=x .1 = 2 J 


5 Chief, Problem Definition and Analysis Branch. 
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Define the terms 


p(E) = PjE' 1 

k+i 

r — i 

or(E) = -h ) (Piaj + pj)r 

k+i 

t(E) = -h 2 Pi } ajE 0 " 

j = 2 

and the characteristic equation (in terms of E) for the system of difference 
equations results from setting the determinant of the matrix 

[P] s [p(E)[l] + ff(E)[A] + t(E)[A3 2 3 

equal to zero. As is well known (see any text on linear algebra) both the 
determinant and the eigenvalues of [P] are identical to those for [P*] where 

[p*] = [b][p ][b ] _1 

for all nonsingular [B]. Furthermore, whether or not [A] has multiple eigen- 
values, [B] can always be chosen so that 

[BJUHB]" 1 = [T] 

where [T] is upper triangular and the elements of its diagonal are, of 
course, the eigenvalues of [A]. But since 

[b][a] 2 [b] = [b][a][b3 -1 [b][a][b3” 1 « [T] 2 

this choice of [B] also makes [P*3 upper triangular. It follows that 

|p| = |P*| = n(p(E) + Aic(E) + A?t(E)) 

where the Aq are the distinct eigenvalues of [A]. This proves the stated 
hypothesis . 


Discussion 

The results just presented show that a predictor -corrector process 
applied to a set of simultaneous, linear, differential equations with constant 
coefficients automatically "detects" the eigenvalues of the differential 
equations and the success or failure of the numerical method is measured by 

1. Its accuracy in resolving the eigenvalue for which it is most 
inaccurate 

2k 


2. Its stability with respect to the eigenvalue for which it is most 
unstable 

Therefore, to study and compare numerical methods as they apply to the solu- 
tion of systems of differential equations (ll) , it is sufficient to study and 
compare them as they apply to a single differential equation. In the case of 
nonmultiple roots, this single equation is simply 

u* = Au + f(x) (35) 

where A may be complex and represents the "worst" eigenvalue of the system. 

That it is generally impossible to estimate the magnitude of the eigen- 
values by the magnitude of the individual elements of a matrix is shown by 
the following examples in which all three matrices have (except for the limi- 
tation of 8 significant figures) the same eigenvalues, -1, -10, and -100. 


-2451-57 52 

9523-1044 

-2225.7133 

\ 

- 765.51692 

2970 • 5032 

- 695.60020 


- 690 . 98884 

2671.4410 

- 629.92802 _ 


”- 41.774674 

65.261307 

15.813105 ~ 


22. 518472 

-39.291311 

-11.264427 


_ 40.012615 

-75-965121 

-29.934014 _ 


”-95.492090 

-H 8.40717 

6 . 8224432 " 


-8.6303246 

-12 . 144652 

0.26226723 


-85.039375 

-114.33994 

-3.3632580 



This can be an important consideration in some programs that attempt to con- 
trol step size automatically using norms based on the elements in individual 
rows or columns (see, e.g., ref. 8). 

Equation (35) is, basically, the representative equation used through- 
out this paper. (For practical reasons, a more convenient expression will 
actually be analyzed, see eq. (37)0 In most papers the representative 
equation is presented as (la), of which equation (35) is a very special form. 
Nearly always, however, commitments about the nonlinear equation are based 
on local linearization for the simple reason that the nonlinear form is 
intractable. In this approach, the reader may regard equation (35) as char- 
acterizing equation (la) when A symbolizes the average value of 8 f/8u 
over some interval (ref. k, p. 207 ), or the Lipshitz constant over the same 
interval (ref. 5, p. 2l6) . 

Accuracy . - For the purpose of studying the accuracy of a difference- 
differential approximation we choose equation (35) as the representative form 
and permit A to be complex throughout the analysis. Equation (35) is still 
unnecessarily general, however, since, for testing the accuracy of a numerical 
method, we can replace f(x) in some interval with its equivalent Fourier 
series and draw out from the latter the term Age^n-* which represents the 
highest frequency in f(x) we wish to resolve by the numerical computation in 
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the chosen interval. All lower frequencies are automatically more accurately 
determined. Hence, for studying accuracy, equation (35) may be replaced by 
the simpler form 


which has the general solution 


u 


1 = Au + Ae^ x 


Ax Ae 
u =■ ce + 


|iX 


p - A 


(37) 

(38) 


where c is an arbitrary constant and both A and \i can be complex. A mea- 
sure of accuracy of any numerical method is given by its worst approximation 
to either term in equation (38) in the presence of the other. 

Stability . - Equation (35) is also used as the representative form for 
analyzing the stability of differ ence -different ial approximations. For such 
studies the term f(x) can be omitted since it does not affect the stability. 
Whether or not a method is stable depends upon the magnitude of the roots to 
the characteristic equation for the difference equations which the method 
generates when combined with the representative form. There is a large and 
important subset to equations (ll) , associated with positive definite forms, 
for which all the eigenvalues are real and positive, that is, for which the 
A in equation (35) is real and negative. For this subset the stability 
criterion can be developed from the simple normalized equation 


u' 


-u 


(39) 


in which the independent variable is real. This form is often used (see 
refs. 2, 3j 9) y and for a numerical method to be stable when applied to equa- 
tions (ll), it must certainly be stable for equation (39). But it is not 
sufficient, and in a later section some methods are shown to be stable for 
equation (39) but not for equations (ll) . 

When the eigenvalues of equations (ll) are complex, the A in equa- 
tion (35) becomes complex. The generalization of the numerical stability 
criterion, under these conditions is presented in reference 10. For an anal- 
ysis that applies this generalized criterion to several numerical methods, 
see reference 11. 


Two ways of approaching the study of the stability problem for complex 
eigenvalues are suggested. One is to consider the simultaneous equations 


w£ = 


w r 


^2 ~ -W^l + 2 v ^2 


(40) 


where v is always real . The advantage of this way is a given set of 
differ ence -different ial equations can actually be used numerically to check 
the results. The necessary and sufficient condition that a numerical method 
be stable for equations (ll) is that it be stable for equations (40) in which 
v is real. When |v| > 1, equations (40) have real eigenvalues. When 
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| v | < 1, equations (40) have two conjugate complex eigenvalues, and the case 
v = 0 represents the condition for pure imaginary eigenvalues. 

A second way to study the general stability problem is to consider the 
normalized equation 

u 1 = e^u (4l) 

which is a special form of equation (35) where 

A =• e iW (42) 

and f(x) = 0. The study of equation (4l) avoids much of the algebra required 
to analyze equations (40) but necessitates the use of complex arithmetic. 
However, in modern computer languages the latter is not a serious disadvan- 
tage. The necessary and sufficient condition that a numerical method be 
stable for equations (ll) is that it be stable for equation (42) for all 
values of u>. 


THE GENERAL ANALYSIS OF INCOMPLETE, MULTISTEP, PREDICTOR, 
ONE -CORRECTOR METHODS 


General Discussion 

The following study applies to all numerical methods that integrate 
ordinary differential equations making use of a predictor followed by just 
one corrector, with both the function and its derivative in the same final 
family. (A study of complete methods is presented in a later section.) These 
incomplete methods include, for example, Hamming r s method, the Adams - 
Bashforth -Moulton methods, Milne 1 s method, etc. In fact, the result of using 
any predictor in table l(a) followed by any corrector in table l(b) can read- 
ily be determined both as to accuracy and stability. The analysis permits 
one to calculate exactly (except for roundoff error) what a digital computer 
will produce after any number of steps when any one of the methods is applied 
to equation (37)* 

T he operational form . - Consider a predictor -corrector sequence forming 
the first family u^ 1 ^ 

k+i 

( a j u n+J + ajhu^+j) , J - k + 1 - j (43) 

j = 2 

and the final family u 

k+i 

u n+k ~ Pi^+k + Pi hu n+k + ^ (Pj u n+J + Pj^ u n+j) (kk) 

Notice that by multiplying the top equation by and subtracting the 

result from the second equation, we can form a third equation 
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,y k+1 

u n+k = Pihu n +k + ^ (Pj u n+J + Pj hu n+j) 


(45) 


j= 2 


This indicates the lack of uniqueness in the values for the coefficients in 
the difference -differential equations for combined predictor -corrector 
methods . 

The following is a simple example of this lack of uniqueness. An Euler 
predictor followed by a modified Euler corrector is generally written in the 
form 

(l) . u ' 1 

u n+i = u-n + hu n 

n 1 (46) 

/ 

u n+i = ^n + | (4+1 + 4) 


Clear ly, the analytical result of such a method is not altered by the modifi- 
cation 


(l) , I 

u n +i - u n + ^ u n 


u n+i 


= c cC u 


(i) 


n+i - u n “ hu ' 


A h r W a f 
n ) + u n + 2 V Un+1 u n/ J 


(47) 


where c D is an arbitrary constant. As an example, setting c G = l/2 gives 
the combination ^ 

C 1 ) ^ , » 

u n+i = u n + hu n 

1 ( ) ■ ( 48 ) 

u n+i = \ ( u n+i + u n + hu n+i ) 

If the corrector is used only once , equations (46) and (48) are identical in 
accuracy and stability when applied to equations (ll) . 

Returning to our development, and choosing, without loss of generality, 
equations (43) and (45), we introduce the operator E and the representative 
equation (37) • This combination produces tne linear difference equations 
which can be written in matrix form 


(49) 


- 

k+i 


— 




k+i 


E k 

-I 

( a j 

+ Aha4)E^ 

J 


u (l) 

u n 


V 

L 

1 J(j.h 
a j e 


j - 2 





= Ahe 1 * 11 " 

j = 2 



k+i 






k+i 


-AhplE k 

i 


+ AhPj)E J 


u n 


I 



j- 2 






A 1 
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Solving for u n , one can divide out from the left column with the result 


k+i 


k+x 


0=2 


(Iaj + AhL 2j + = hAe^ ^ (Rij + AhR 2j )e J ^ h (50) 

J 


where 


L lj 

= Pj 

Ln 

L 2 j 

= ajPi + Pj 


L 3j 

= atp’ 

R 2j 


- 


(51a) 


Equation (50) is the operational form of multistep, two-iteration methods and 
the L and R are the coefficients in the operational form. Coefficients in 
the operational forms of a variety of methods are given in table II. 


Equation (51a) can be inverted. That is, the coefficients in the 
difference -differential equations can be expressed explicitly in terms of the 
coefficients in the operational form, provided Rn 0. Thus 


a j - ( 1*2 j “ R lj)/ R ll 

a j = Raj/Rn 
Pj = L ij 
Pj = 


(51b) 


Fortunately, cases for which = 0 appear to have little practical use. 

This inversion is the key to the construction of optimum numerical methods , 
since both stability and accuracy depend fundamentally only on the operational 
form (discussed later) , not on the difference-differential equations. We now 
seek only to analyze given methods, not to develop new ones. 


Before proceeding, it is convenient to introduce the following two 
definitions . 


DE(E) = the coefficient of u n in any operational form 

IJthn 


NU - (p - A) times the coefficient of Ae H 
operational form 


in any 


(52) 


In particular, for multistep, two-iteration methods 


DE(E) = E* 1 


k+i 

^ (Lij + AhL 2 j + A 2 h 2 L 3 j)El J 
0=2 


(53a) 
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k+i 


(53b) 


NU = h((i 


A) ^ (Rxj + AhR 2j )e Jtlh 


i=i 


The particular solution- - Referring to the section on the operational 
solution of difference equations, we can immediately write the particular 
solution in the form 


k+i 


V (Eu + AhE ZJ )e J|ih 


1-1 


u 


nr- 


kp.h 


k+i 

^ (L X j + AhL 2 j + A 2 h 2 L 3 j ) e J|lh 




or , using the definitions in equations (53 ) , 


% = 


MJ 

^ ~ A / DE(e M ' h ) 


(5*0 


The complementary soluti on.- The explicit evaluation of the complementary 
solution requires a knowledge of the roots to the characteristic equation. 

In the case under consideration, the characteristic equation is 


E k 


k+i 

^ (la j + AhL 2 j + A 2 h 2 L 3 j)E J = 0 


or, using equation (53a) , 


DE(E) = 0 o (55a) 

Let the roots to equation (55) be such that 

(E - Ai)(E - A 2 )(E - A 3 ) - * • =■ 0 (55h) 

The complementary solution is then 

u Uc = Ci(Ai) + g 2 (A 2 ) + 03(7^3) + • • • (56) 

where Ci>c 2 ,C 3 , . . . are constants fixed by the initial conditions, and, 
conventionally, Ai is the principal root while A 2 ,A 3 , • • • are all spurious 
roots introduced by the particular numerical method. 


30 



Accuracy 


The error in the particular solution . - Comparing equations ( 38 ) and (5^)> 
one can now derive the error in the particular solution. Defining the error 
term 


Exact particular solution of difference equation 

— - - - ^ i =— = • - - * * — 1 

Exact particular solution of differential equation 
it follows that 


( 57 ) 


— - MJ - DEfe^k) 

^ " DE(e^ h ) 

Introducing the values of HU and DE(e^' 1 ) into the numerator of (58), and 
collecting coefficients of Ah, gives 


( 58 ) 


k+i 


~k+i ~ 

W 1 ** * £ (Li 3 + 

+Ah 

^ (Lgj - R ld + phR 2j )e Jiah 

j =1 


- 


er n = 


DE(e M ' h ) 


( 59 ) 


The coefficient of the A 2 h 2 term is zero since R 2 j = L 3 j (see eq. (51a)). 

Equation ( 59 ) is the exact error which any incomplete multistep, two- 
iteration numerical method makes in calculating the particular solution to 
equation ( 37 )* However, since all methods being studied here are polynomial 
approximations, it is more significant as well as more convenient, to express 
the errors in powers of h. 


Setting 


qh 


(H h )' 


z! 


one can show 


Z-o 
k+i 


k+i / k+i k+i ' 

DE(e^ h ) = 1 - Iaj + h/jik - n ^ Jlaj - A ^ L 2 j 
j-s \ <U 2 j ~ 2 ) 


+ 0(h 2 ) 


(60) 


As will be shown presently (see eq. (72)), to have any polynomial fit, at all 
the equalities 
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Lij - 1 


k+i 

I 

0=2 

k+i 


yi [l 2 j + jLi j ] - 


0=2 

must hold. Under these conditions 

k+i 


DE(e^ h ) _ 


k+i 


lim 
h-»o h 

Next, noting 

one can show 
where 

er, 


= (u - A) ) (o - DLij = (u - A) 


J 20 


0=2 


0=2 


(uh) Z+1 J Z V z ^ h) J 
n L v. 

1=0 1=0 


Z T Z- 1 


^Fjj. ~ 


er m + er M.g _ er(i 
h(h - A) h(p - A) 

■ k+i 


Hi 


(uh) 

Z ! 


z | X [zjZ " lRl j + jZli J ] " k *| 
- k+1 

^ (0 - l)Iaj 


0=2 


er, 


Ah(qh) 


(i2 


^ [ZJ Z_1 R £ j + J Z (L 2j - Rij)]j- 


k+i 

£ (3 - DWj 


0=2 


(61) 


(62) 


(63) 


m 


( 65 ) 


Now it is important to notice that erj^ is a "global” error; that is, 
it is the precise error in the particular solution at a given x no matter 
how many steps are taken. The local error, or error in the embedded poly- 
nomial, such as erp in equation (4), is repeated n times after n steps. 
Because x = nh, the polynomial error can grow, for a given x, to 
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approximately (x)(er p )/h. Hence, rather than er^, the term erq represents 
the accuracy of the local polynomial in a given method . 7 


If a method is to represent a polynomial approximation of order L to 
the particular solution between n + k - 1 and n + k, then the terms inside 
{} in equation (64) must sum to zero for all l =• 0, 1, 2 , . . . , L and the 
terms inside {} in equation ( 65 ) must sum to zero for all 
l = 0, 1, 2, . . . , L - 1. This gives at once the equations that must be 
satisfied by the coefficients in the operational form for a method to provide 
a polynomial approximation of a given degree to the particular solution. 
Specifically, 
k+i 


I 

j =i 


[z(k + 1 - o) Z_1 E 


1 3 


(k + 1 - jflaj] = k 1 


* = 0 , 1 , 2 , 


L ( 66 ) 


k+i 

^ [ Z (k + 1 - j) Z_ 1 R 2 j + (k + 1 - j) Z (L 2j - Rxj)] = 0 , 

j =1 Z = 0, 1, 2, . . L - 1 (67) 


We shall see later (in the discussion of the complementary solution) that the 
fulfillment of these conditions guarantees a polynomial fit of degree L to 
both the complementary and particular solutions. 

Notice that equations ( 66 ) and ( 67 ) are independent of A and q so the 
coefficients of R and L can be tabulated once and for all. Table III does 
just that for any one- through five -step, predictor, one -corrector method. 

The table can be used both to provide the conditions that a given polynomial 
be embedded in a method, and, with equations (64) and ( 65 ) , to find the error 
in the result. Similar tables were used to calculate the error er^ for the 
methods presented in table II. 


The error in the complementary solution . - Just how well the complemen- 
tary solution is approximated by a numerical method depends upon how close 
the principal root, A x , lies to e^ 1 , since the analytic solution is propor- 
tional to e^ = ( e ^h) n , and the numerical solution is proportional to (Ai) . 
Let er-^ be defined by 

er^ ~ Ai - e^ 1 (68) 

Substituting e ^ 1 in equations (55b) and (55a) 

k+i 3 

JAh V 

e L 

j=g 

7 The term (q - A) in equation ( 63 ) can be misleading. Recall that it 
comes from the expansion of DE(e |J '' :1 ) in powers of h. As q -» A, the term 
0(h 2 ) in equation ( 60 ) takes over, and the error does increase. When q = A, 
the solution of the representative equation degenerates and the analysis 
proceeds along different lines . 



er^(e ?sh - A 2 )(e Ah 


- A3) 


= 
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or 


er A “ 


-DE(e Ah ) 

.n (e Ah - Aq) 

1-2 


(69) 


Again we are 
k 

value of II (e^ 1 
i=2 v 


interested in finding only the lowest 
- Aq) for h = 0 can be determined as 


order error term, 
follows : 


The 


since the principal root must equal one at h = 0. This, in turn, reduces to 





Note the similarity with equation (62) which appears in the denominator of 
er^. Putting equations ( 52 ) and (70) in equation (69), and expanding in 
powers of Ah, one can show 


er A 



+ ZJ Z_ 1 L 2 j + l(l - l)J z " 2 L 3 j] - k 


Z -2t 


k+i 


^ (j - DLij 


j=2 


(71) 


The first nonvanishing term in (71) determines the order of the error in the 
complementary solution. (Notice this is a "local" error, comparable to the 
term er^. In fact, one can show that er^ = er^ when n = A.) 

If a method is to provide a polynomial approximation of order L to the 
complementary solution between n + k - 1 and n + k, the terms inside {} in 
equation (71) must sum to zero for all l = 0, 1, 2, . . ., L. This derives 
the equations 


3 ^ 



r 


k+i 

^ [(k+l-j) Z L 1<3 +l(k+l-j) Z_1 L 2 j+z(Z-l)(k+l-j) 2 " 2 L3j] = k Z 

l = 0, 1, 2, . . L (72) 

Equations ( 72 ) are independent of 7\ and \i so, again, the coefficients of L 
can be tabulated once and for all (for one- through five -step methods, see 
table IV), By means of such tables and equation (7l)> the errors in the com- 
plementary solution of predictor, one -corrector methods can be readily 
determined. Examples are presented in table II. 

Discussion of accuracy .- One can show (by using R 2 j = I< 3 j) that equa- 
tions ( 66 ), ( 67 ), and ( 72 ) are not independent. Hence, as was mentioned 
earlier, the satisfaction of equations ( 66 ) and ( 67 ) is sufficient to guaran- 
tee a local polynomial fit of degree L to both the particular and comple- 
mentary solutions. Thus two sets of conditions must be fulfilled to assure a 
given accuracy for a predictor, one-corrector process. If another corrector 
is added, another set of conditions must be met, etc., as will be shown later. 

On the other hand, if we wish to present the conditions for a predictor 
only (or for a corrector only which has somehow been brought into balance) , 
we can set the aj and c xj terms in equations (43) and ( 51 a) equal to zero. 

(in this degenerate case the (3 values become the coefficients of a predictor 
with the condition 0i = 0.) Then equation ( 67 ) is an identity, and Lij-Pj, 
Rij = Pj. Equations ( 66 ) and ( 72 ) are identical, and they both amount to 
successive columns in the table preceding equation (4) adding up to zero 
(or, in equation (4), er p = 0 for l = 0, 1 , 2, . . . , L) . In this case the 
first unmatched term for erp, x or er^ is the same as the "error constant" 
term given by Henrici on page 223 in reference 5* 

In summary, equations ( 66 ) and ( 67 ) are the most general conditions for 
accuracy imposed upon the coefficients of two difference-differential equa- 
tions forming a predictor -corrector sequence under the conditions: 

(a) The differential equations are of the form given by equations (ll). 

(b) Polynomial approximation is used. 

(c) The difference-differential equations are of a form represented by 
equations (43) and (45)* 


Stability 

A general discussion of stability is given in the last section of the 
report where the Dahlquist criterion is extended to cover all linear, numer- 
ical, quadrature formulas with combined effects of predictors, correctors, 
modifiers, etc. In this section the specific stability of predictor, one- 
corrector methods is considered. 
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Some definitions , - To begin with, let us examine some terminology com- 
monly used in discussions of numerical stability. A set of ordinary differ- 
ential equations is inherently unstable if the real part of one or more of the 
eigenvalues of equations (ll) is positive. There are two classes to consider. 

Class 1. The initial or boundary conditions are such that the analytical 
solution grows exponentially. Then the numerical solution must also 
grow exponentially and is, therefore, by definition, unstable — but it 
is not necessarily inaccurate. If these inherently unstable differential 
equations are transformed into difference equations, the latter are rela- 
tively inaccurate (sometimes referred to as "relatively unstable," see 
the discussion below on induced instability) if any of the spurious roots 
are greater in absolute value than the largest principal root- 

Class 2. The initial or boundary conditions are such that the destabi- 
lizing eigenvectors are suppressed and the exact analytic solution has 
no terms which grow exponentially. Under these conditions the numerical 
solution will still eventually increase exponentially, usually because 
of small truncation errors that excite the inherent unstable terms, but, 
if for no other reason, because of errors brought about by roundoff. 

This can be a particularly vicious form of instability and its control 
requires methods outside the scope of this paper. 


When a set of ordinary differential equations is reduced to a set of 
ordinary difference equations, the latter have an induced instability if the 
real parts of all the eigenvalues of the differential equations are all less 
than or equal to zero, but one or more of the eigenvalues of the difference 
equations has an absolute value greater than one . This instability is 
obviously associated with the particular form of the difference-differential 
equation chosen for the computations and can, therefore, be- controlled. 


Remembering that the eigenvalues of the difference equations are func- 
tions of h, we find two ways of providing this control. One is to develop 
methods that are stable when h = 0. Then for small enough 8 h, the method 
is always stable. The Dahlquist stability theorem applies to this study. 

The other aspect is to develop methods that are stable for as large a value 
of h as possible. In order to discuss this problem, we will refer to a 
stability boundary, | Ah [ c ^ which defines a critical value of h, any increase 
of which causes the absolute value of one or more eigenvalues in the differ- 
ence equation to exceed unity. Thus a method has induced instability when 

A = e iw q 


Tt 


/2 < w < ir 


( 73 ) 


Ah > Ah c \ 


avoid argument, we either omit the possibility of neutral stability 
at h = 0, or further require that all eigenvalues that lie on the unit 
circle when h = 0 move into it as h starts to increase. A method having 
the latter property is defined by equations (120) . 
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A discussion of the significance of this stability boundary is presented in 
a later section entitled "The relationship between accuracy and stability." 

A final remark concerns the case when the real parts of all the eigen- 
values from the differential equation are less than zero and the absolute 
values of all the eigenvalues from the difference equations are less than one, 
but some of the difference eigenvalues are greater in magnitude than the dif- 
ferential eigenvalues. This condition is sometimes referred to as being 
"relatively unstable." If we are to maintain a consistent definition of sta- 
bility, this terminology is misleading. Such cases are stable since all solu- 
tions, as well as their errors, approach zero with increasing x. They may 
possibly be relatively inaccurate, however, the possibility occurring when the 
f(x) terms in equations (ll) are negligible and the offending eigenvalues are 
the least heavily damped. In any case one should be cautious about trusting 
the numerical solution of the asymptotic decay of a function when the level 
falls beneath the product of the truncation error of the method and its 
maximum resolved amplitude. 

The unit circle . - Since the induced stability of a method depends upon 
the magnitude of the roots to the characteristic equation DE(E) = 0, a quick, 
visual, representation of the stability of a method is displayed if we plot 
all the roots to the characteristic equation on a complex plane. For the 
construction of these plots let 0 < h < 1 and 


where 0 < u> < it . Describe a unit circle with center at the origin of this 
plane. In the range 0 < to < it/2 "the differential equation is inherently 
unstable and the principal root, at least, must fall outside the circle for 
h > 0. In the range it/2 < to < jt the differential equation is inherently 
stable and in this range any point falling outside the circle presents a value 
of h and A (in complex form) for which the numerical method has induced 
instability. 

The two extremes: Adams -Moulton methods , Milne methods . - The accuracy 

of most of the popularly used predictor -corrector formulas has been compro- 
mised to avoid induced instability. That this compromise can be resolved in 
tnar$ ways is evident from the number of predictor -corrector methods in com- 
mon use. However, all of these formulas lie between two extremes which we 
will refer to as the Milne methods and the Adams -Moulton methods. These 
extremes have certain identifying features which appear immediately in the 
stability plots described above, but which can also be traced to coefficients 
in the operational form and even to coefficients in the difference- 
differential equations themselves. 

The Adams -Moulton methods are defined as those for which all the spurious 
roots fall on the origin when h = 0. The principal root at h=0 must 
always, of course, fall on the intersection of the unit circle and the posi- 
tive real axis. When h = 0, the characteristic equation (see ( 53 )) reduces 
to 
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k+i 

E k - l 1j eJ = 0 

j=2 


The necessary conditions that make the spurious roots zero and the principal 
root one are 

Li 2 = 1 

bij = 0 , J - 3j • • • , k + 1 

In such a case the characteristic equation reduces to just 

(E - 1 )E k_1 = 0 

which clearly satisfies the conditions. Inspecting equations (51b), we see 
this means 


3a = 1 

3j “ 0 , j ~ 3 ? • • • > k + 1 

Thus, in Adams -Moulton methods, only the value of u nearest the corrected 
point is given a nonzero weight in the corrector. 

In this report methods which have at least one spurious root on the unit 
circle when h = 0 are referred to as Milne methods, and those which have all 
the spurious roots on the unit circle are referred to as total Milne methods . 9 
In the latter case the characteristic equation becomes 

k+i 

E k - Y L i t j E ' J = ( E - 1) .6 (E - e 10 3) = 0 


and the are determined such that the L x j are real and the accuracy is 

optimum. In these cases Lq ^ 0 and, from equations (51), 3k+i ^ 0* 

Thus, in the total Milne methods the value of u farthest from the corrected 
point is given a nonzero weight in the corrector. 

All multistep methods, with induced stability at h = 0, lie between 
these two extremes. On the one hand they have no spurious roots on the unit 
circle so they are more likely to be stable for nonzero h than the Milne 
methods. On the other hand, if some of the roots at h = 0 are permitted to 
lie off the origin - but still within the unit circle - they gain some free- 
dom which can be used to choose the L x j, j - 2, 3, . . .,k + 1 so that they 
will be more accurate than Adams -Moulton methods with equivalent step number. 


^hat are called total Milne methods in this report are referred to as 
optimal methods in reference 5* The author prefers the former description 
since, in practice, stability troubles generally prevent such methods from 
being optimum. 
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The relationsh ip between accuracy and stability # - To appreciate the role 
that induced instability plays in assessing the merits of various methods, we 
should distinguish carefully between the requirements for stability and accu- 
racy. The error in any method, as it is given by the lowest order, nonvanish- 
ing, truncation terms (e.g., the values for er^ and erA in table II), loses 
its significance when these terms exceed about one tenth. In fact, to rely 
on such error estimates, the step size should be chosen so that 

(||ih|, | Ah| ) max <0.1 

in methods applied to the integration of 

u T = Au + Ae^ X 


Now one usual requirement for predictor -corrector methods, programmed 
for general use, is that the induced stability boundary |Ah| c , see (73)^ be 
as large as possible, generally greater than 0.6. The question immediately 
arises: Of what value is a method that is stable in a region where it is 
inaccurate? The answer is supplied if we consider the application of the 
method to the integration of simultaneous equations. Let us consider a 
predictor -corrector process for which 1 7\h | c = 0.6l. To prohibit induced 
instability and give an accuracy measured by the values of er^ and er^ per- 
taining to the method, it is certainly sufficient that 


( | \ij_h j , | cr^hj ) <0.1 , for all i (7^a) 

where the \i ^ and cr^ are determined by the differential equations (see 
eqs. (ll) and (l6)). Although (7^-a) is a sufficient condition, it is not a 
necessary one. The necessary conditions for both stability and accuracy (in 
the sense used above) are that 


(1) ( | M-ih | , |cTih|) max < 0.1 , for all i representing those q and a 
one seeks to calculate to the specified accuracy 

( 2 ) l a i h lmax < l Ah lc > for a11 i 

( 7^b) 

This distinction between accuracy and stability is sometimes important. 10 
Perhaps the easiest way to describe the situation for those unfamiliar with it 
is to give an example. Consider the equations 

= - 1 . 38 w! - 0.8lw 2 


with the initial values 


Wg “ “2.l6wj_ - 1*92^2 

w^o) = -2.9905 

w 2 (o) = 4.0010 


(75a) 


(75b) 


10 For example, the so-called "stiff" equations arising in the study of 
nonequilibrium fluid flow. 
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For many purposes the step size h - 0.2 is perfectly acceptable for solving 
this problem in the range, say, 0 < x < 10, using a predictor -corrector method 
for which |Ah| c = 0.6l. For the equations as presented, the reason is cer- 
tainly not obvious. But if the eigenvectors of the equations are formed by 
the relations 


Wl = 5U! - 3 u 2 
w 2 = 10u! + 4 u 2 

one can show equations (75) are equivalent to 

ui = -3^1 

/ 

u 2 = —0 . 3 u 2 

Ui(0) = 0.0001 

u 2 (0) = 1.0000 

The analytic solutions of these are, of course, 

u x = 0.0001e -3x 


with the initial conditions 


Up — l.OOOOe 


— O • 3X 


(76) 


(77a) 


(77b) 


Now suppose we are not interested in values of u x and u 2 when they are 
< 0.0001. As we have seen ; the stability of predictor -corrector methods, when 
applied to simultaneous equations, depends upon the "worst" eigenvalue of the 
differential system, in this case - 3 . Since -(-3)(0.2) = 0.6 < 0.6l, stabil- 
ity is assured. This corresponds to the second condition in (74b) . As far 
as the required accuracy is concerned, stability is all that is necessary for 
Ui, since it is smaller than the allowable error to begin with and, being 
stable, cannot grow. The accuracy of u 2 will be acceptable for any method 
that makes | er^J < 0.0001 for |Ah| = 0.06. This corresponds to the first 
condition in (74b) . 

Of course, if the initial conditions for equations (75a) are changed to 

wi(0) = 2.0000 

w 2 (0) = 6.0000 

the situation is quite different. Now the analytic solutions are 

Ux = l.OOOOe " 3X 
u 2 = l.OOOOe "°‘ 3X 

and in the interval 0 < x < 3*1* Ui > 0.0001. To get the same resolution as 
before, the step size would have to start at 0.02. By the time x ^3.1, how- 
ever, a step size of 0.2 would again be satisfactory. 



Clearly, the importance of these considerations depends upon the size 
and complexity of the set of differential equations being studied. We can 
calculate, however, the most one can gain by using the necessary conditions 
(74b) rather than the sufficient conditions (74a). Let |?\a| represent the 
largest eigenvalue or frequency (i.e., | A a | =• (|a|, \[i \ 1 ^) we wish to 
resolve in a given problem. Let the maximum 11 step size be given by 
|hA a | “ O.O 3 . Consider the two curves in sketch (c). They represent the 

maximum step size for a given A 
determined on the basis of accuracy 
( 0 . 03 / | A | ) and stability ( | Ah | c / | A | ) , 
where | Ah | c is the stability bound- 
ary of the numerical method. If A c 
is the negative eigenvalue largest in 
magnitude and we wish to resolve it, 
then | A a | > |A C | and the accuracy 
curve governs the step size, giving 

h min # Howev * er J if I A al < l A cl it 
may be possible to increase the step 
size until it is governed by the 
stability curve, giving h,^. It 
is impossible to increase h further 
by the methods discussed in this 
report. Under these conditions the 
ratio of the maximum to the minimum 
step size is 

Me 

0.03 

which has a significant effect in the solution of lengthy problems. 

Some examples .- Next, we construct the unit circle in the complex plane 
and plot in the same plane the roots to the characteristic equations for a 
variety of predictor -corrector methods. In each example the roots are calcu- 
lated for a step size equal to zero and the locus of these points is indicated 
by flagged symbols. The step size is then incremented by 0.1 and each root is 
recorded accordingly. The number of points plotted varies because of scale 
limitations, but in no case does | Ah | exceed one. The value of u) in the 
representative equation u T = e 1Cx) u is varied through the range jt > to > jt /2 
to study induced instability, and through 0 > to > - at/2 to study the region 
of inherent instability. 

The location of the principal root under these conditions is shown in 
figure 1. This root must, of course, be common to all methods, and a measure 
of the accuracy of a method is displayed by how closely one of its set of 
points falls on those in figure 1 . 

1:L Simple two-step methods with er^ ^ er^ 0.0l(hA a ) 4 are available. 

If |hA a | = 0.03 such methods have an error 8l(l0) “ 8 . Smaller step sizes 
in these (or more accurate) methods in a machine using eight -place floating 
arithmetic would be useless. 


(78) 



Sketch (c) 


X max 


hi 


min 


4l 




Figure 2 represents the method resulting from the combination of a four- 
step, Adams -Bashforth predictor (row 5j table l(a) ) followed by a four -step,, 
Adams -Moulton corrector (row 4, table 1(b)). All of the spurious roots lie 
on the origin at h = 0, forming a triple root there and causing the variation 
of the roots with h to be quite different for small h (SAp/dh -*■ ch -2/3 as 
h -*■ 0, i f l) from what it is^for h w 0.1 and higher. The method presents no 
stability problem until 1 7\h| >0.7. Beyond this value a spurious root exceeds 
unity for A = ±i. hot ice that the upper left-hand circle (A = -1) shows the 
method is stable for Ah even less than -1.0 for real eigenvalues. This is a 
case for which a stability analysis that makes use of only real A would give 
quite erroneous results regarding the value of a method for general 
simultaneous equations. 

Figure 3 displays the root structure of Hamming's method without modi- 
fiers (row 6, table l(a), and row 6, table l(b)). Four roots are involved 
and two of the three spurious roots do not fall on the origin at h = 0. When 
A = -1, one of the spurious roots starts at 0.422 and proceeds in a positive 
direction along the real axis for increasing h. It crosses the principal 
root at h « 0.265 causing a degenerate instability there. For Ah < -0.5 it 
falls outside the unit circle and the method becomes definitely unstable. 

This is a case for which the general stability boundary is determined by 
considering only real eigenvalues. 

The method Hamming finally proposed, and the one usually programmed and 
referred to by his name, uses two "modifiers . " An analysis and discussion of 
modifiers is given in the next section, where it is shown that Hamming's 
method with modifiers is equivalent to predicting with row J, table l(a) and 
correcting with row 7, table l(b) . The roots to the characteristic equation 
for this method are shown in figure 4. The root structure is completely dif- 
ferent from that shown for the unmodified method. There are now five roots. 
One of the four spurious roots lies on the origin when h = 0. The_ eigen- 
value limiting the stability is now complex, occurring when A = e 2 ^/ 3 and 
the stability boundary is | Ah| <0.6, slightly greater than the unmodified 
one. 

Typical of what can be done to increase the stability of these four- and 
five -root methods is shown in figure 5. This figure illustrates the roots to 
the characteristic equation for a method proposed by Crane and KLopfenstein 
(ref. 11) which amounts to using the four -step predictor in row 8, table X(a) , 
and the four -step Adams -Moulton corrector (row 4, table 1(b)). For this 
method, all of the spurious roots fall in the unit circle for |Ah| <0.9. 

The classical Milne method (row 6, table l(a) , and row 5j table 1(b)) has 
the characteristic roots shown in figure 6. Four roots are generated; two of 
the spurious roots lie on the origin and one at -1 when h = 0. (So, in our 
terminology, it is a Milne method, but not a total Milne method.) It is well 
known that the corrector alone is unstable for all negative Ah. Chase 
(ref. 3) showed that for real Ah the combined, predictor -corrector process 
is stable for -0.3 > Ah > -0.8, a conclusion which is represented here in the 
upper circle on the left. Clearly, however, a glance at the entire left 
column in the figure shows that for almost all complex eigenv al ues with 
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Figure 3 «“ Hamming's method without modifiers 



Figure k.- Hamming's method with modifiers. 


oO°co a 







Figure 5.- Method of Crane and Klopfenstein, 




negative real parts, the Milne predictor -corrector method is unstable for all 
h > 0. This result is also presented in reference 11. 


The necessity for considering the stability problem in the entire complex 
plane has certainly been demonstrated. We will present one more example, 
however, since it provides a background for some of our subsequent discussion. 
Consider the total Milne method composed of the two-step predictor in row 10, 
table I (a), followed by the Milne, two-step corrector (row table 1(b)). 

This method has been proposed by Stetter (ref. 9). Since now both the predic- 
tor and the corrector have only two steps, only two roots appear in the char- 
acteristic equation. They are shown in figure 7> the left and right columns 
of the previous figures being collapsed into two circles for convenience. 

For all real Ah the method is stable for 0 > 7\h > - 1. It is the most 
accurate conventional (i.e., incomplete) predictor, one-corrector method 
that can be devised for arbitrarily small h (er^ and er^ are given in 
table II (k))^ and may be used if one is sure that all the eigenvalues of the 
differential equation are real. As the left circle in figure 7 shows, how- 
ever, like the classical Milne method, it is unstable for all imaginary 
eigenvalues . 



METHODS WITH MODIFIERS OR NOimJHDAMEHTAL FAMILIES 


Incomplete, Predictor, One -Corrector Methods 
With Modifiers 

Only fundamental families were used to construct the methods studied in 
the preceding section. It is quite possible, of course, to hold in memory, 
and subsequently use, combinations of u and u* which were calculated in a 
previous cycle of computation but are not members of the final families. The 
equations relating these combinations are often referred to as modifiers. In 
the terminology of this report, using a modifier corresponds to constructing 
a family that is not fundamental. The principal purpose of this section is 
to show (by operational methods) that modifiers are artificial in the 
following sense: 
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Any method 12 with modifiers or nonfundamental families can be 
identified with a method without modifiers composed only of 
fundamental families . 


One of the simplest types of modifiers is that which weighs only past 
families of the dependent variable and not its derivative. It requires no 
further iterations, but it does, of course, require more memory. The set of 
difference-differential equations with modifiers that we will now analyze can 
be written 

(i) 


k+i 


u n+k = Y, ( a 0 Un+ J + ha J u n+j) 

(i) 

n+k 

, ( 2 ) ' \~ ' . 1 t , 

hPi^n+k + > (3j u n+J + hPju n+ j) 


j=2 

u n+k = T 2 u n+i + t 3 uA+ 
k+i 


( 2 ) _ _ ..( 1 ) , _ „( 3 ^._ 1 


, ( 1 ) 

+ T 4 U n +lc-l 


(3) 

u n+k 




_ - ( 3 ) ( 1 ) (3) ( 1 ) 

u n+k ~ a i u n+k + cr2 u n+k + ^3^n+k-i + tf4 u n+k-i 


( 79 a) 

( 79 b) 

(79c) 

( 79 a) 


Equations (79a) and (79 c ) are the conventional predictor -corrector equations 
studied in the previous section. Equations (79b) and (79<3-) are the modifiers, 
weighing previously calculated, nonfundamental families by the constants tj 
and o j • Applying these equations to the representative equation and intro- 
ducing the operator E, one derives the matrix equality 


E* 


-T 2 E k - T 4 E k 1 


_k -nk-l 

-cr 2 E - a 4 E 


0 


E 


-AhPiE k 


-t 3 E' 


E 


,k-i 


-aiE k - cr^ -1 


where 


k+i 


-C 

0 

-Cf 

E 

k+i 


a 


,k 


' 

r u (*>‘ 

u n 


Fa 


( 2 ) 

u n 

= Ahe^ hn 

0 


( 3 ) 




L u n _ 


0 


( 80 ) 


C a - Y (aj + Ahaj)E d 


Fa = 


a je 


* J\ih 


k+i 


j=2 

k+i 


Cp = 


(Pi + MiPi)E J 


Fr = 


Q « Jph 

pje 


. 1=2 


j =1 


(81) 


12 We only consider incomplete methods but the same conclusions apply to 
complete ones . 
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If we solve equations (80) for u n , we can divide out E k_1 from the first 
and third columns, and E^ from the second column. There results 


_k+2 


k+i 

-I 

j=o 


k+i 


E 


Lmj (Ah) 


m-i 


‘ u n = hAe 


phn 


Jph 


Rl(Ah) 


m-i 


m=i 


i = -l 


m=i 


( 82 ) 


where L* and R* are fairly simple combinations of the constants cc, t, P, 
and a in equations (79)* It is clear that, if we set 


equation ( 82 ) 


k e = k + 2 
* 

Lm,j+2 = Lmj 
R m, j +2 = Rmj 

J e =k e +l-j 

can be written 


(83) 


kp+i 


E 


j= 2 


kp+i 


E 


L mj (Ah) 


m-i 


u n = hAe 


phn 


m=i 


= JePh 


Rmj (Ah) 


m-i 


(84) 


. 1=1 


m=x 


Except for the subscript e, this is the same as equation (50),, the opera- 
tional form of the difference -differential equations (43) and (45) , composed 
only of fundamental families. Thus the k-step method with modifiers given 
by equations ( 79 ) can be identified exactly with a higher step method without 
modifiers. 

More general forms of modifiers, weighing more past families of the func- 
tion as well as its derivative, can be analyzed. They would simply increase 
the number of terms with powers of E in the square matrix in equation (80) 
and lead to characteristic equations in (82) of higher order. By substitu- 
tions similar to those in ( 83 ) , the final operational form can be again iden- 
tified with equation (50). This correspondence of methods (as they apply to 
linear differential equations) is always established when the operational 
forms are identical. 


Hamming 1 s Method With Modifiers 

A good example of equivalent methods, one using two modifiers, and the 
other using no modifiers but having one more step is given by analyzing 
Hamming 1 s method as it is usually programmed (see ref. 2). Hamming 1 s 
modified method can be written 
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( 1 ) 

u n+4 

(2) 

u n+4 


= u n + I h (zu-n+3 - u n +2 + 2un+i) 


u 


( 3 ) = 

■n+4 


C 1 ) 112 r (1) ( 3 K 

u n+4 - 12I V*n +3 " Un+3 J 
1 r r (2) ' 1 

g |_9 u n+3 - u n+i + 3h\^ u n+4 + 2u n+ 3 


- u n+2^)_ 


u, 


L n+4 


121 


r (1) 

|_9u n +4 + 


( 3 )~ 

112u n+4 


( 85 ) 


By a straightforward calculation, using the formulae in the previous section, 
one can show that these have the coefficients in an operational form given by 
table Il(c). Using equations (51b), we immediately find two difference- 
differential equations that have the same operational form. These are 

u n+s = u n+4 + Un+i - Un + j h(2u n+4 - 3 u n+3 + 3u n +2 - 2u n+ i) 
u n +s = |^ 126 u n+4 - lUu n+2 + 9 u n+i + h(^42u n+5 + 108 un-H 4 “ 54u n+ 3 + 24u n + 2 ^j 

( 86 ) 


Equations (86) represent a conventional, five-step, incomplete, predictor - 
corrector method composed of two fundamental families which, except for round- 
off errors, gives results identical to those obtained using Hamming 1 s modified 
method when applied to equations (ll) . 


Discussion 

Consideration of the previous sections raises the question as to the 
nature of the relationship between families and steps. In what was just pre- 
sented, a five-step method was duplicated by a four -step method with addi- 
tional families. How far can this be carried? The answer is that any k-step 
method can be reduced to a one-step method if the number of families is 
increased appropriately. (The converse is not true; the minimum number of 
families is the number of iterations in incomplete methods and one plus the 
number of iterations in complete ones.) The next question that arises is 
whether or not this introduction of families serves any really useful purpose. 
After all, any given method basically evaluates a polynomial of a certain 
degree using an amount of data stored in memory necessary to attain that 
degree. From this point of view, there is little difference between a method 
expressed with modifiers and the same method reduced to fundamental families. 
Possibly, a few storage locations can be saved and a few arithmetic or logic 
manipulations eliminated by using one or the other. Most likely, roundoff 
accumulations will differ, but these are not considered here. 
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There is quite another point of view, however, from which the family 
concept can play a valuable role. One can show that, if the proper families 
are constructed, any polynomial method can be reduced to a several -family , 
one -step method, the step size of which can be changed at will after each 
advance. To derive such constructions systematically, however, requires a 
theory that falls outside this report. 


COMPLETE MULTISTEP PREDICTOR -CORRECTOR METHODS 


Introduction 


We define complete predictor -corrector methods as those in which the 
final values of the function and its derivative are not members of the same 
family. A k-step, two-family, complete method is represented by the two 
equations 


(i) _ 

u n+k “ 


u n+k 


ft ». C 1 )' 

Pihu n+k 


k+i 


C a j u n+J + 


J=2 

k+i 


Y C p j un+j 


- (l) 

Pju n+ J 


. 1=2 


a jhu n +j J 
Pjhun+J 


(87) 


This method requires only one iteration per step. Thus, if the calculation 
of the derivative dominates the computing time, when compared to incomplete 
methods, complete ones can 

(1) Use twice as many steps for the same amount of computing time, or 

(2) Cover the same interval with the same step size in half the 
computing time. 


We seek to find whether or not, on the basis of accuracy and stability, these 
gains are real or fictitious. 


Notice 

this report 
retained in 
be 


that equations ( 87 ) are not composed of (what has been defined in 

to be) fundamental families, since both u n+ j and u n+ j are 
memory. Complete methods that use only fundamental families would 


( 1 ) 


u n+k 


k+i 

Y C a J Un+J + a 5 hu n+J ) 

j=2 


u n+k 


, (i>* 

= Pihu n+ k 


k+i 

+ y C p 3 un+j + 

j=2 


pjhun+J ) 


( 88 ) 
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for one iteration, and 



(2) 

u n+k - 



t ( 2 ) 1 

u n+k = 7i u n+k + 7 ihu n+k 


k+i 

+ ^ C 7 j Un+ J + 7 J hu n+J ) 

j=2 


for two iterations- We show in the next section that almost all operational 
forms given by equations (87) can be constructed from some form of equa- 
tions (88) . 


Analysis and Discussion 


Introduce the representative equation (37) into equations (87) and there 
follows the matrix equation 


k+i 

k+i 




k+i 

E k - ^ (aj + Ahctj ) 

- Z 


A 1] 


Z 





= Ahe^ htl 

j=2 

k+i 

k+i 




k+i 

-AhpiE k - Y ( F j + Ahp j ) E J 

K*- £ 


u n 



j=2 

j=2 




j= 2 


(90) 

Compare equations (90) and (49) and the difference between the incomplete and 
the complete methods begins to appear. The left column in (49) contains the 
term E^ which can be factored out, leaving a characteristic polynomial of 
order k. In equation (90) no such factoring can be made, the left column 
is completely (hence the terminology) filled with terms E°, 

and the characteristic polynomial is now of order 2k. 


Solving equation (90) for u n gives the operational form 



2 k+i 


2k+i 


£ (Iaj + AhL 2 j)E 2 * + 1 -jL n - hAe^ hn £ 


1=2 


1=1 


(9D 
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where L and R are combinations of the a and p. In the simplest case when 
k = 1 these combinations are 


1*12 ” a 2 + 02 
L13 = a 2 + a 2Pi 


^22 “ p2 a 2 “ a 2p2 


Rll “ Pi 


^12 “ P2 “ a 2Pl 

Rl3 - p2 a 2 “ a 2p2 


Ij23 ~ P2 a 2 - a 2p2 


Although there are exactly seven terms on both sides of the equations, invert- 
ing them, so as to express a, p in terms of L,R, would be difficult because 
of their nonlinear form. 

If we use only fundamental families, equations (87) reduce to equa- 
tions (88), and the combinations 


cla + PiQH + \ (Pi a i+i-i " a iP- 


l 3 + vPi^j+i-i - ^iPj+i-iJ 


“ Pi 9 ^lj “ Pj > J 2 , 


are formed. But these equations can be inverted, in general, since 

Pi = R n > Pj = R ij > Pj “ L ij > j = 2 , . . ., k + 1 ( 94 a) 


a j + ^ii a j + ) (Hi jdj+i -i - Lija^) - L 2 j , j - 1 , . . . , 2 k + 1 


the latter being a set of linear simultaneous equations for a. For example, 
if k = 2, equations ( 94 ) reduce to 


Pi = R11 > P2 “ R i2 j P2 - L 12 


P3 = R13 9 P3 " L13 
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(95b) 


R 11 

0 

Rl2 

Hu 

Rl3 

Rl2 

0 

Rl3 


1 0 

-Li 2 1 

-Lx3 -Li 2 

0 -L13 




L 22 

a 3 


1*23 

t 

CX»2 

=. 

1*24 

t 

a 3 


I»25 


respectively- In this inverted form L and R are arbitrary and, in partic- 
ular, L24 and L 2 5 can be equated to zero. This leads to a set of equations 
from which a and P (coefficients in the difference-differential equations) 
can at once be derived for any combination of L and R on the left side of 
equations (93 ) > provided only that the determinant of the matrix in equa- 
tion (95b) is not zero. In other words, except for the restriction just given, 
any operational form contained in equations (87) for k = 1 can be identified 
with an operational form from equations (88) for k = 2. This situation 
remains true for higher values of k, so that equations (88) are sufficiently 
general to represent almost 13 all complete, two-family, predictor -corrector 
methods. Further, the relations between the coefficients in its operational 
form and the coefficients in the difference -differential equation can always 
be inverted by means of equations (94), provided only that the determinant of 
(94b) is not zero. 

The accuracy and stability of the two-family, complete methods can be 
studied using the same analysis as that presented for equation (50). For the 
complete case, simply set I©j and R 2 j equal to zero. It is apparent from 
equations (50) and (91) that a k-step complete method will have accuracy and 
stability features associated with a 2k-step, incomplete method. Superfi- 
cially, this appears to violate the Dahlquist criterion. When the latter is 
applied to operational forms, however, we see that in reality no such viola- 
tion occurs. The difference in stability between complete and incomplete 
forms is discussed later. 


Examples 


An example of a complete, two-step, predictor -corrector method contained 
in equations (88) is 


%l +2 = ^ [>*1 + l8u n + h(l5un+x + 7u^ ^)] 


h /- (l) ' n (l) ' (l) ' 

u n+2 _ u n+i + J 2 \5 u n+2 + 8u n+1 - u n 


(96) 


13 Cases in which equation (94b) are overdetermined have not been 
investigated. 
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The coefficients in its operational form are 




L mj 



R mj 

m 

2 

3 

4 

" 5 
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4 
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1 

24 
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0 

0 

10 

6 

-2 

0 

0 

2 

55 

-59 

37 

-9 







Divide by 24 


er H = 1^2 - A)H 4 , er A ~ 0(H S ) 

and the error terms, referenced to the computing step h, are 

er pi = ik > er P 2 = d 3 Ah 4 , er A ~ 0 (h 5 ) 


(97) 


Equations ( 96 ) employ a two-step predictor and an Adams -Moulton, two- 
step corrector; and they have Adams -Moulton type stability. An incomplete, 
two-step, predictor -corrector method using the Adams -Moulton corrector can be 
written 


(1) 

u n +2 


= -2 . 8u n+1 + 3-8u n + h( 3 • 4i%+i + l-4u n ) 


u n +2 - u n+i + 3 ^ ( 5n n+2 + 8 un +1 


(i)' 


u n 


(98) 


and the coefficients in its operational form are 


\ j 

Lmj 

Rmj 
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3 

1 

2 

3 

1 

12 

0 

5 

8 

-1 

2 

-6 

18 

0 

17 

7 

3 

17 
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Divide by 

12 



er p = - A)H 4 , er x ~ 0(H 5 ) 
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The lowest order error terras for the method given in (98) are identical to 
those given in equations (97) when referenced to the computing step h. (To 
order h 4 this is the minimum possible error if an Adams -Moulton corrector 
is used.) We see that equations (96) have the same (lowest order) accuracy 
as equations (98) , but require one less it eration per step . 


For a fair comparison, the error terms in all methods should be refer- 
enced to the representative step size H, defined as the distance the solution 
is advanced by two iterations. For equations (96) we have 


and for equations (98) 


h = | H 
h = H 


A remark is in order here with regard to the adjustment of er^ and 
er^ to conform with the representative step. Since we are now concerned 
with the efficiency of a method as it applies to the complete calculation of 
the differential equations, the global, rather than the local, error should 
be used as a measure. This will account for the fact that if one method uses 
half the step size of another, it commits its local error twice as many times. 
Hence, if in general, h = Ha 


and 


er |j.(k n ) -* ^ er 1 j i (a n H n ) 
er A (h n ) i er^(a n H n ) 


(99) 


Errors referenced to H are given for the two methods in tables follow- 
ing equations (96) and (98)* Clear ly, on the basis of accuracy, the complete 
method is to be preferred, having l/8 the error of the incomplete method with 
the same corrector. Since both have Adams -Moulton type stability, all the 
spurious roots vanish at h = 0. There renj&ins the question, however, regard- 
ing the magnitudes of the spurious roots for h ^ 0. The characteristic 
equations DE(E) = 0 for the two methods are 

24E 4 - (24 + 55Ah)E 3 + Ah(59E 2 - 37E - 9) =0 (100) 

for the complete method and 

12E 2 - (12 - 6A h + 17A 2 h 2 )E - Ah(l8 + 7Ah) = 0 (101) 

for the incomplete one. We at once see the complete method has two more 

spurious roots than its counterpart if h ^ 0. The magnitudes of the roots 

for real and complex Ah are shown in figures 8(a) and 8(b) . In terms of h, 

the calculation step size, the complete method has induced instability when 
| Ah | > 0. 3> whereas the incomplete method has no induced instability until 
| Ah | = 0-5- However, if we again refer our measurements to the representative 
step size, the boundaries are |AH| = 0.6 and |Ah| = 0-5, respectively. 
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(b) Incomplete, two-step, two -iteration method given by equations (98) . 
Figure 8.- Concluded. 


In summary, if the evaluation of the derivative dominates the computing 
time, then, for the same computing time, a complete method with the two-step, 
Adams -Moulton corrector is more stable than, and has about l/8 the error of, 
an incomplete method with the same error. The statement regarding error is 
based on only the first term in a series expansion. Experience indicates it 
is not reliable for | Ah| >0.1. 

GENERAL ANALYSIS OF INCOMPLETE MULTISTEP METHODS 
WITH MULTIPLE CORRECTORS 


Derivation of the General Solution for a Fixed Corrector 

The process defined by equations ( 43 ) and ( 44 ) can be generalized such 
that, during the same cycle of computation, an arbitrary number of correc- 
tors -- that is, an arbitrary number of iterations -- is used. Even when 
each of the correctors is different, the general solution can be derived by 
the technique outlined below. The solution is quite complicated, however, 
and does not appear to be of practical interest. The special case, when two 
correctors, both different, are used, is analyzed later in this section. If 
all the correctors are the same, the general solution has a rather simple 
form; and it provides us with the ability to study the effect of the number 
of iterations on the accuracy and stability of a method. 

The equations for the fundamental families of an incomplete method in 
which m - 1 correctors are used, all identical, can be written 
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If 

FT 1 — 

II 

k+i 

1 

( a j u n+J + ha j u n+j) 


0= 2 


( 2 ) , (l)' 

u n+ k = hPxUn+k 

k+i 



(Pj u n+J + hPjUn+j) 


<1=2 



( 102 ) 


» (m-i) 

u n+k - kPi u n+k 


k+i 


(0j u n+J + h Pj u n+j) 


j=2 


Using the notation defined in equations (8l) we introduce the representative 
equation and the operator E and derive the matrix equation 


E k 0 0 . . 
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0 -b^ E k 
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= Ah^ hn 

F P 

-bxE 11 

E k - Cp 


_ u n 


F P_ 


b x = Ahpi 


(103) 


(10k) 


Expand the determinant about the right -hand” column and the characteristic 
equation simplifies to 

T>T1 / T - lN t_,1c « , m-i _ r , m-s . m-3 1 

DE(E) — E — Cccb^ - Cp[b^ + b^ +. . .+1] 

_k 


- E- - Cab?" 1 - Cp, VkF 1 


= 0 


(105) 


Introduce the notation 

L 1 = (cy + Ahaj)bx _1 + (P i + Ahpj) 


, m-i 

-n t -i m— i 0 r 1 - bn 

Ri = a^jbi + P-s * k - 


1 - b 


m-i 


1 - bn 


1 - b X 


(106) 
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and the operational form turns out to be 

k+i ^ k+i 

E^ - 

Z_. u I " t-, 

j=2 J j= s 

If 


Ahe^ hn ) R^* 1 


k+i 


E k - ^ L j E<I = ■ A l)( E - A 2) * • * 


k+i 


J=2 

the complete solution can be written 


u n = CiCAi ) 11 + C 2 (A 2 ) n + . . . + Ahe ^ 1 * 111 ^ 


R J e 


Jnh 


k+i 

kuh V T Juh 
e - > Lje" 


(107) 


(108) 


The case m = 1 represents one iteration per cycle of computation* that 
is* a predictor without a corrector. This follows from equations (106) since 
Pj and Pj disappear from Lj and Rj when m = 1. On the other hand* if 
rn -+ oo* the equations are independent of a,j and aj (provided | b x | <1* which 
is a necessary condition for the convergence of the iterations), and the 
corrector equation in its implicit form emerges. 


A Discussion of Some Simple Predictor -Corrector Methods 


If we use equation (108) to inspect the simple predictor -corrector scheme 
(an Euler predictor followed by a modified Euler corrector) 

(x) t 

Un+i = u n + hu n 


X s ) = 


h f (x) 


u n+x = u n + 2 Un+1 


+ u. 


n 




h 


u n+x = u n + 2 l u n+i + u n 


we find the complementary solution to be after m iterations 

,m+3T ,n 


u nQ _ Ci 


1 + 


Ah _ 


i 

x 2 


(109) 
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This result serves as an excellent example of the danger of analyzing a 
corrector, ignoring the effect of a predictor. If the modified Euler method 
is studied alone -- as an implicit identity — one finds the complementary 
solution n 


Un C 


= Ci 


" l + (Ah/2) 
-1 - (Ah/2)- 


( 110 ) 


instead of equation (109) • This shows at once that the modified Euler method 
is stable for all negative values of Ah (see ref. 12). But clearly, as m 
becomes large, equation ( 109 ) reduces to equation (110) only when jAh| <2. 
Hence, when used as a corrector in a predictor -corrector sequence , the modi- 
fied Euler method"Xor trapezoidal rule as it is sometimes referred to) is 
violently unstable for Ah « -2. 


A further study reveals that er^ and er^ for the Euler -modified-Euler 
method behave in the following manner for increasing numbers of iterations in 
a cycle of computation 

m erA 

1 -h 2 (q 2 /2) -h 2 (A 2 / 2) 


2 -h 3 (n 3 - 37v- 2 )/12 -h 3 (A 3 /6) 


3 h 3 (|_i 3 /l2) h 3 (A 3 /l2) 


h 3 (n 3 /l2) 


h S (A 3 / 12) 


This method is often used in programming the "method of characteristics" in 
the study of hyperbolic partial differential 14 equations. It is of interest 
to notice: 


1. The order of error in the predictor is one less than the order of 
error in the corrector, but the order of error in the method is the 
same as that for the corrector after one application of the 
corrector. 

2. The coefficient of the h 3 term is improved by a second application 
of the corrector. 

3. The error, as measured by the lowest order in the truncation, is not 
further reduced if the iterations are continued. 

If the Euler predictor is replaced by the Nystrom predictor (row 2 in 
table 1(a)), the error sequence with iterations is 


l4 When more than one independent variable is involved, the reference 
step H should be redefined. 
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m gtjj. er A Spurious root 

1 -h 3 (|j. 3 / 6) -h 3 (A 3 /6) -1 + Ah 

2 h 3 (|J. 3 /l2) h 3 (A 3 /l2) --|Ah 

oo ii 3 (n s /l2) h 3 (A 3 /l2) Ah 

Here, we see that the predictor by itself has the same order of error as the 
corrector but is less accurate, and, further, is unstable. One application of 
the corrector (a total of 2 iterations) 

1. Increases the accuracy such that no improvement on the coefficient 
of h 3 is made by continuing the iterations 

2. Produces a stable numerical method. 

Basically, this process is the one used in several computer programs to solve 
for the flow in front of blunt bodies travelling at high speeds. 


Incomplete Multistep Predictor Two-Corrector Methods 

Next, the incomplete methods previously analyzed are extended by adding 
one more corrector with arbitrary coefficients. Only a brief sketch of the 
procedure is given here, mostly for the sake of thoroughness, since the prac- 
ticality of using two correctors is open to question. However, the added cor- 
rector has a decidedly stabilizing effect. In fact, it is shown that, for 
two-step incomplete cases, a stable, two-corrector method can be constructed 
having error terms one order higher than is possible for stable, one -corrector 
methods. 



where, by the argument used to derive equation (45) Pi can be made zero with- 
out loss of generality; but, by the same argument, the term containing jx 
cannot. Although a weighted value of the middle equation (ill) added to the 
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final equation can be made to cancel the term involving (similar to what 

was done in equation (44)), a term containing Un+1 would then appear 
amounting to a method different from that given by equations (ill) with 
7 1 = 0 . 

Paralleling the one-corrector case studied in a previous section, one can 
construct a matrix equation and derive the operational form. There results 

r k+i 4 ^ k+i 3 

E k - )T Ljnj (Ah) m_ 1 (-u n = Ahe ^ hn ^ e J[lh ^ R^CAh )" 1 ' 1 (112) 

i i=2 m=i I j=a m=i 


from which 


and 


k+i 


DE(E) = E k - y E J L m j(Ah) 


ra-i 


(113) 


J-s 

k+x 


m=i 


WU — h(|x - A) ^ e ^ ^ R m j(Ah) 

= 0 ) 


m-i 


(114) 


j =1 


m=i 


For j - 2 , 

3 , . • . 

(since Ln, L 2 i, 

. . . 


Li j 

= Pj7x + 




L 2j 

= a-jPbx 

+ 0J7 1 + 

Pj7i + 


L 3j 

= a jPx7i 

+ a jPx7x 

_ 1 f 

+ Pj7x 


L 4j 

t t r 

- a t jPx7x 



and for j = 1 , 

2 , • - • 






■ + 




R 2j 

= a jPx7x 




R 3 j 

= L 4j 



These equations 

uniquely determine 

L and R, 

the co 


(115) 


tional form, for given a, 0 , and 7 . Once again, the accuracy of the method 
is represented by equations ( 63 ) and ( 68 ), and the stability depends on the 
magnitude of the roots to DE(E) = 0. 


Equations (ll5) can be inverted if y± is set equal to zero. Thus 
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(116) 


a j - ( L 3j - R 2 j)/k 2 i 
0j = (L 2 j - Rij)/Rn 
7J = L lj 


ai =■ 


M 

Pj 

7j 



and difference-differential equations that represent any operational form for 
which En and R 12 are not zero can at once he written. The practical 
consequence of these limitations is not known. 


To find the error in the particular solution, the numerator of (58) is 
expanded in powers of Ah. In this case R3j = L 4 j and the coefficient of 
the A 3 h 3 term vanishes identically. Choosing the L and E so the coeffi- 
cients of the (Ah)°, (Ah) 1 , and (Ah) 2 terms are zero gives, respectively, 


k+i 

^ [l(k + 1 
j =1 


j ) 1 + (k + 1 ~ j) Z Lij] - 'k 2 = 0 , 

l - 0, 1, 2, 


L (117a) 


k+1 

^ [ Z(k + 1 

j =1 


+ (k + 1 - j) 2 (L 2j - Rij)] = 0 , 

Z = 0, 1, 2, 


L - 1 (117b) 


k+i 

I 

-1-X 


[ Z(k + 1 - j) 2 Raj + (k + 1 - j) 2 (L 3 j - Rgj) ] = 0 , 

1 = 0 , 1 , 2 , 


L - 2 ( 117 c) 


which are the accuracy conditions for equations (ill) . One can show that, as 
in the single corrector case, both the particular and complementary solutions 
are fit locally by polynomials of order L if the above conditions are 
satisfied. 


The leading terms in the series expansion for er|j, are determined by 
evaluating for Z = L + 1 the remainders in the expressions -- they are no 
longer equalities — (117) multiplying these remainders by (ph) lrl ’ 1 /(L + 1) t , 
Ah(ph)^/L'. , and A 2 h 2 (ph)R“ 1 /(L - 1)1, respectively, and dividing each by 
k+i 

y (i - l)Iaj as in the development of equations (64) and (65). 
j* 2 


The derivation of the error in the complementary term is identical to 
the derivation of equation (71) ■ Thus one can show 
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I * k+1 


Jj YT-< Y, + + Z ( Z Z(Z-l)(Z-2)j^L4j]-k 2 


ex^- 




k+i 

£ (o-l)Li^ 




(118) 


Finally, to compare fairly with other methods, all error terms should be 
expressed in terms of the representative step H where 


h = 3H/2 

and the adjustments are made according to equations (99) . 
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Example . - The incomplete two-step, double-corrector method given by 
1 


M 

u n+2 

U (2) 

u n+2 


u n+2 


= ^ [ -248u n+1 + 337% + h(302% +1 + 124%) ] 


1 

375 


-432u n+1 + 807u n + h ( 89%+2 + ?88un+i + 305% 


h f ( 2 ) * - t A 

“ Un + 3 (^ u ti+2 + 4u n+1 + u n ) 


( 120 ) 


has the operational form represented by 


V 

L jm 

Rjm 

2 

3 

1 

2 

3 

1 

0 

1125 

375 

1500 

375 

2 

1068 

1182 

89 

788 

305 

3 

540 

642 

0 

302 

124 

4 

302 

124 





Divide by 1125 


6% = (0.028p 2 - 0.04QJ.A - 0,020A 2 )(Jl 3 H 5 
er A = -0.032A 5 H 5 


The accuracy of the method as measured by the computational step size is given 


,e% = (0.0057P 2 - 0.0080V - 0.0040A s )p 3 h 5 
er^ - -0.0063A 5 h 5 


( 121 ) 


and as measured by the representative step size as shown in the table. 
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There are two roots to the characteristic equation 


1125E 2 - E(l068Ah + 540A 2 h 2 + 302A 3 h 3 ) -1125 - ll82Ah - 642A 2 h 2 - 124A 3 h 3 = 0 

( 122 ) 

and they are shown in figure 9- The spurious root starts on the unit circle 


Negative A Positive K 







Figure 9 .- Stability of predictor, two-corrector method given by equations (l20) . 


at -1. Although it is impossible to tell from the figure, a close examina- 
tion of the results shows that the spurious root moves into the unit circle 
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f° r all complex A given by A = e 1C0 for which 7r/2<u)<3jt/2. Further, 
the spurious root remains inside the unit circle for |Ah| < O.57, or 
| AH| < O.38 which establishes the stability boundary according to equa- 
tion ( 73 )* This serves as an example of a total Milne method (all the spurious 
roots are on the unit circle when h = 0) that has no induced stabilities for 
a range of h > 0 . The stability boundary can be increased to at least 
| Ah| c = O.56 if the maximum error is allowed to increase by a factor of about 
3 / 2 . 


COMBINED RUNGE-KOTTA AND PREDICTOR -CORRECTOR METHODS 


Introduction 

Now let us consider the overall result of increasing the number of iter- 
ations in a cycle of computation. In the analysis of the equation u ! = F(x,u) , 
it is clear that each successive iteration requires the re-evaluation of the 
term F(x,u) at a value of u different from any used in previous iterations. 
As the number of iterations increases, appropriate choices of a, J 3 , 7, etc., 
permit us to match the final result with a Taylor series expansion of F(x,u) 
in the u direction through any given order -- regardless of step number. In 
other words, the accuracy of fit to the complementary solution depends on 
both the step number and the number of iterations; and its series expansion 
can be matched with arbitrary accuracy by increasing either the one or the 
other independently . 

Next consider the particular solution to the differential equations. 
Conventional predictor -corrector formulas are constructed using a fixed and 
equal spacing of the independent variable, x, a spacing we have designated as 
h. In such cases the number of samplings of F(x,u) in the x direction is 
determined entirely by the number of steps used in the method. No new infor- 
mation regarding the variation of F(x,u) with x is supplied by increasing 
the number of correctors in a cycle of computation. Since, in general, the 
error associated with a method must depend on its worst fit to either the 
particular or the complementary solution, equispaced, predictor -corrector 
methods are limited in accuracy by their step number, regardless of the number 
of iterations. 

Clearly, if both u and x are varied in the successive correctors, the 
complete series expansion of F(x,u) in both x and u can be matched to an 
indefinite order, and arbitrary accuracies to both the complementary and par- 
ticular solutions obtained. The development of this concept leads, at once, 
to equations that merge the predictor -corrector and the Runge-Kutta methods. 


On the General Form of the Equations 

One can write a set of difference -differential equations that represent 
a complete, combined method with k steps and m iterations. The result 
would be a set of formulas in which all families of all the values of the 
function and its derivative calculated in a cycle of computation and held in 
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memory are weighted with appropriate values of a , p, 7 , etc. (see, e.g., 
ref. 13) . By considering only the two-step case, we will see the complexity 
involved in such a completely general expression. For the practical purpose 
of actually constructing methods with fixed accuracy and stability, much 
simpler subsets of these general expressions appear to be satisfactory. The 
real crux of the problem of constructing optimum methods depends upon whether 
or not a set of equations can be found that explicitly relates the coeffi- 
cients in the difference-differential equations to the coefficients in the 
operational form (e.g., eqs . (51b), (94), and (ll 6 )). This question, in turn, 
depends upon how the operator E is brought into the matrix equation from 
which the operational form is constructed. 

The discussion of combined methods is more readily presented if we 
consider simple sketches in which the following symbols are used: 

□ predicted value of the function, u ^ 

A derivative of the predicted value, du^ 1 ) /dx 

O corrected value of the function at a point previously predicted, u 
<^> derivative of the corrected value, du/dx 

On the definition of step size in combined methods. - Using the symbols 
defined above, let us construct a cycle of computation for a two-step, 
predictor -corrector method (e.g., eg. (98)), and the standard, one-step, 
fourth-order, Runge-Kutta method (eqs. (136)). 

When presented graphically, as in the sketches, it may at first appear 
that singling out a length h and calling it step size is a rather arbitrary 
procedure. In fact, an ambiguity in the use of the words "step size" and 
"step number" can easily arise when predictor -corrector methods are combined 
with Runge-Kutta methods; although, as we shall presently see, the term can 
be given a unique and quite natural definition that applies to the two dif- 
ferent approaches, either individually or in combination. 

In conventional predictor -corrector schemes, such as that shown in 
sketch (d), the function and its derivative are calculated at equispaced 
points only, and a value of the function is (or can be without loss of accu- 
racy) "output" at each point. The spacing is quite naturally referred to as 
the step size and the resulting step number can be and is used as the funda- 
mental parameter in describing both the accuracy and stability of the method. 
See, for example, the Dahlquist stability theorem. On the other hand, in 
Runge-Kutta methods the function and its derivative are calculated at points 
other than those at which it is most accurately represented or intended to 
be output. In sketch (e) it is shown at the midpoint, but this is totally 
unnecessary. Nevertheless, regardless of the number of intermediate points 
used, or their spacing, the Runge-Kutta methods are always referred to as 
one-step methods. 
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This definition is in keeping with the usage throughout this report and in all 
references about numerical methods of which the author is aware, including the 
recent ones on combined methods (refs. 13 , l6) . 

With such a definition, step size is still a useful parameter for describ- 
ing the accuracy of a method, but it_ no longer has any basic meaning with 
regard to stability . In short, stable combined methods exist which embed 
polynomials of arbitrarily high order regardless of step number . This is 
further discussed in the last section of the report entitled Stability. 

When comparing the accuracy and stability of combined methods, then, 
either with themselves or other methods, one must do so on the basis of con- 
trolled values of such things as : 

1. The amount of memory required 

2. The amount of arithmetic required 

Examples of a two-step method . - In the following, we present a variety 
of two-step methods subjected to the following restraints: 

1. A memory of --at most -- four values of the function and/or its 
der ivative, 

2. The calculation of -- at most -- two families in a cycle of 
computation. 

With the addition of one more family or word of memory, the accuracy and 
stability of any method can probably always be improved. 


Method I. Two iterations, incomplete, uncombined 
tl) _ 

U n+2- 


a 2 U n+l + a 3 U n +h(c, 2 U n+l + a 3 U n> 


u n+2 = £' hu n+2 +/3 2 u n+ l + /3 3 u n + h ( ^2 u n+l + ^3 u 't. > 


Maximum order 
of error for 
stable method 


Data used in corrector 
is encircled 


New values calculated 
in a cycle of 
computation 
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Method 2a Two iterations, incomplete, combined 



Method 2b Two iterations, incomplete, combined 



n 


Method 3 One iteration, complete, uncombined 

U n+2 = a 2 U n + l + a 3 U n +h ^ a 2 U n+l + a 3 U n ' ^ 

u n+2~ “n+2 + #2 u n + 1 + ^3 u n + h </ 3 2 u n + 1 + ^3 u n ^ 


Method 4a. One iteration, complete, combined 



Method 4b. One iteration, complete, combined 

(I) - ^ (l) , . , / (l)' , / (i/ , 

u n+r~ a 2 u n+l + a 3 u n+r-| + h( a 2 U n+r-| +a 3 U n+ r - 2 ) 


u n+2 ' 


=^[ui l / r +^ 2 u n+l +^ rH + h(4u>7'- l +^iV ) r '- 2 ) 



Meth od 5. One iteration, incomplete, uncombined 

u n+2 = a 2 u n + l + a 3 u n +h(a 2 u n+l + a 3 u n ) 


er=0(h 3 ) 



Notice that, in each case, a cycle of computation is completed when the 
integration has been advanced a distance h, and the number of iterations 
refers to the total number of evaluations of the derivative in this cycle- 
At first glance, it appears that the two iteration methods do not belong in 
the same group with the one -iteration methods because more arithmetic is cer- 
tainly required to evaluate the derivative of the corrected function and this 
violates condition 2 of this section. However, this effect is taken care 
of by referring the stability and accuracy terms in all methods to the refer- 
ence step size H- With this important qualification, all methods 1 through 
4 can be compared on the basis of very nearly equal logical complexity. 
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storage and computation time., provided only that the calculations required to 
evaluate F(x,u) in the equation u f => F(x,u) are equal to or greater than 

those required to find u^+r an & u n+2 in a cycle of computation. 


Method 1 represents the classical predictor -corrector sequence; for 
example, an Adams -Bashforth predictor followed by an Adams -Moulton corrector. 
One can show by applying the analysis presented in a previous section that no 
method of this type is stable for arbitrarily small but nonzero h and for 
arbitrary complex eigenvalues if the local error is 0(h 5 ) . Method 2a is 
typical of the combined methods proposed in references l4, 15, and 1 6 . At the 
beginning of a cycle of computation the function and its derivative are cal- 
culated at the point n + r, and the derivative is used to repredict the func- 
tion at n + 2. Neither the function nor its derivative is retained in memory. 
This scheme can be used to develop stable methods with a local error of 0(h 5 ) 
(see eqs. ( 156 )). The order of the error cannot be further increased since 
0 (h 5 ) is the highest order possible for any method with a five term corrector, 
and we are limited to correctors with a maximum of five terms by the condi- 
tions assumed. We can undoubtedly improve the magnitude of the error and 
increase the stability if we make methods (l) and (2) complete. However, this 
would necessitate the addition of another family -- again violating the 
assumed conditions. Method 2b falls into the same class as method 2a, but it 
retains in memory the value of the derivative of the function at the inter- 
mediate point, rather than the value of the function at n. A method of this 
type is studied under equations (159) • It is the most accurate method of all 
those illustrated, having the leading error terms erp, - l/720(pH ) 5 and 
er-^ = l/720(AH) 5 . Its induced stability boundary is \M c ^ 0.3* 

Next, we consider one -iteration methods which can be made both complete 
and combined under the imposed conditions. The simplest one-iteration process 
is the incomplete, uncombined one composed of a single predictor, method 5 in 
our case. The accuracy is severely limited by the requirement of stability. 

If this process is made complete, as in method the accuracy of stable 

methods can be increased by an order of magnitude -- see equations ( 96 ). If 

it is further combined with the Runge-Kutta techniques, as in method 4a, the 
error term for stable methods can be reduced to 0 (h 5 ) just as in method 2 a 

(see eqs. ( 152 )). Another choice of data for a one -iteration, complete, com- 

bined method is presented as method 4b. Rather surprisingly, however, this 
choice is always unstable if the error is to be 0(h 5 ). Proof of the latter 
is given in this part commencing with equations ( 153 ) • 


A Special Class of Multistep, Multi -iteration 
Combined Methods 

Let us next study a class of multi -iteration, combined methods. Consider 
a set of equations formed by using a memory of the function and its derivative 
at only equispace intervals, but, during a cycle of computation, predicting 
and correcting at several arbitrarily placed intermediate points. The anal- 
ysis of such a process is rather simple and shows the connection between the 
classical predictor -corrector methods and methods for multiple numbers of 
steps and iterations. 
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Assume that values of u and u* 
equispaced points hn, h(n + 1), . . 


Memory 



have been computed or given at the k 
. , h(n + k - 1) as in sketch (f) . We 

seek to advance the solution one step 
to the point h(n + k) . We permit 
ourselves up to four iterations, or 
three correctors, but we allow for 
the possibility of evaluating the 
function and its derivative the first 
three times at the points h(n + r x ), 
h(n + r 2 ), and h(n + r 3 ), where the 
r values need not be integers or 
even lie in the interval between 
n + k - 1 and n + k. A set of 
expressions for four fundamental 
families can be written: 


(i) 

u n+ri - 


( 2 ) 


u n+r 2 “ 


k+i 


I 

(aju n+ j + ctjhun+j) 

<5=2 


k+i 


n* U)’ V 

Pi u n+r! + ) 

( 3 j u n+J + Pjhun+j) 

<5=2 



(s) 

Un+r3 


2 1 k+i 

^ 7j4>+rj + ^ (l j u n+J + 7g hu n+j) 

<5=i <5=2 


u n+k 


3 



.. (•))' 
6ju n+rj 


k+i 



( S j%i+J + 


Sjhun+j) 


(124a) 


( 12^+b ) 


(l24c) 


(i24d) 


Equations (124) require four iterations per cycle of computation, 
the following analysis applies directly to three -iteration, incomplete 
combined methods if 


a j ~ a j - 3i ~ 7* - 8* - 0 


} 


All 


and families 2 and 3 are replaced by 1 and 2, respectively. The two- 
iteration case follows by further reduction from the top. 

The incomplete, uncombined, predictor -corrector methods discussed 
previously are obtained from equations (124) by setting rj = k. Classical 
Runge-Kutta methods result if the memory is set equal to zero for j > 2. 
Specifically, the standard, fourth-order, Runge-Kutta method results when 
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(1) k = 1 

(2) r x = r 2 = ~ , r 3 = 1 

A detailed discussion of this method is presented in the next section of this 
part . 

The matrix equation for the combined, incomplete methods defined by ( 124 ) 
follows exactly as it was developed in the previous parts. Extending the 
notation in equations (8l) to include terms with 7 and 5 , we find 


E' 


r l 




-AhpiE' 


-Ah7*E ri 


-AhS^E 1 " 1 


^2 


E 


-Ah72E r2 


-Ah6gE r2 


0 

0 

?3 


-G 


a 


- c 3 


E' 


-Ah&3E r3 


E - C 5 


u. 


( 1 ) 


n 


Ui 


(2) 


■n 


u. 


(s) 


■n 


u. 


n 


= Ah^ htl 


*a 

+ Fp 


ihri 


2 


+ F-. 


j =1 


£ - F & 


j =1 


( 125 )' 


Now if we seek only u n , each term E J common to a column can be factored 
out; and, r egar dless of the choice of r j , the operational form has only 
integer exponents of E. Using straightforward algebraic manipulations, one 
can derive the expressions defined in equations (52) • 


DE(E) = E^ 


k+i 5 

E J ^ L m j (Ah) m-1 
j=2 m=i 


(126a) 


NU = h(|j. - A) 


k+i 


(Ah) 


4 


4-m 

"- 1 V * V e* 1 ) R m .i(*h) 


m-i 


m=x 






The coefficients in the operational form are 

R h - 6 j 

R21 = &37l + &2P1 
R22 = 63/! 

R31 = &37M 


m=i 


(126b) 


( 127 a) 
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(127b) 


R ij = 8 3 

R 2j = ajst + Pj6l + 7jSt 
R 3 j = ^(p?&I+7i^) + Pj7l&3 
R 4J = aj P&I&3 


L ij “ 8 j 

L^j " Hi j + aj5i + Pj&2 + 7j^3 

I»3 j = H 2 J + a j(Pl&2 + 71 ^ 3 ) + Pj72^3 

L4j = H 3j + ajpj 7 g63 


(127c) 


They can be inverted (provided R31J R22 an ^- R13 
expressions 

= R31/R22 

7 x = (BSi - pM 2 )A 

7 I - Rla/R^s 


are not zero) to form the 

(128) 


= R 


* 

ij 


a j = ( L4 j 

Pj = (Lsj 

7 j = (L 2J 


- R 3 j)/Rli 

- R 2 j - a,jRia)/R22 

- Rij - aj-Rli - PjR*2)/R*3 


(129) 


8 j 

- Iaj 


a J 

= L5 j/ R 3 l 


H 

= (Rsj - 

a jR21 ) / R 22 

1 

73 

= ( R 2j “ 

ajRh - PjrL)/ r : 

8 i 




( 130 ) 


Again we use equation (58 ) to calculate the error in the particular 
solution. Expanding the numerator }.n powers of Ah we find the conditions 
for making the coefficients to (Ah) J vanish for j = 0, . . .,3- Since 
Lsj = R 4 j y the coefficient to (Ah) 4 is identically zero. 
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3 


k+i 

^ Z(k+l-j) Z Ri j+(k+l-j) Z Li j +yzr^ Rij-k Z =0 , 
j =2 . 1=1 


k+i 


1=0,1, . . . , L (131) 


y Z(k+l-j) Z " 1 R 2 j+(k+l-j) Z (L 2r R lJ ) ^ rjRjj-ZrJ^Rgj =0, 

0=2 j=i 

l = 0,1, . . . , L - 1 

k+ 1 


(132) 


l (k+l-j ) 1 _1 R 3 j+(k+l- j ) 1 (I© j -R 2 j ) 


-j ^ “- L - r) * . 

r j R 2j-^ R 3 j 


*0, 


0=2 


k+i 


j =1 


l = 0, 1, . • L - 2 (133) 


y Z (k+l-j) 2 _1 R 4 j+(k+l-J ) 1 (L 4J -Raj ) 


- riR 31 =0, 


j=2 


1 = 0 , 1 , . . ., L - 3 ( 13 U) 


These are the accuracy conditions for equations (124) . Once more one can show 
that both the particular and complementary solutions are fit locally by poly- 
nomials of order L if the above conditions are satisfied. Except for the 
addition of one more term (which is multiplied by A 3 h 3 (qh)^” 2 /(L-3) I ) the 
value of er^ is determined just as in the discussion under equation (lrfe). 
The error in the complementary solution is 



<3=2 


(k+l-j) Z L x<j +Z (k+l-j ) 1 _1 L 2 j+Z (Z -l) (k+l-j ) 1 " 2 L 3 j 


+z(z -1) (z -2) (k+l-j) 1 _3 L 4 j+Z(Z -1) (Z -2) (z -3) (k+l-j ) 1 “ 4 Ls j 


-k* 


'k+i 


( j-l)Li J 


0=2 


(135) 


To compare with other methods all error terms should be expressed in terms of 
the representative step H where 


h = 2H 


since four iterations are made to advance one step h. 
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Accuracy of the Standard Fourth -Order Runge-Kutta Method 


The standard Runge-Kutta formula (see, e.g., ref. 17) can be written in 
the predictor -corrector notation used in equation (124) . It is 


u (i) 

u n+o.5 

( 2 ) 

u n+o . 5 


( 3 ) 

u n +i 


u n +i 


- u n + 2 ' 1Un 


1 , ( 1 ) ' 

u n + 2 ^ u n+o .5 


= u, 


n 


+ 4 h 


%i + 3 


h 


(2)' 
n+o . 5 

~ ( 3 ) ' f (2) » (1) ' 

u-n+i + 2(u n+0 .5 + u n+0 .5 


+ u n 


(136) 


If - h is replaced by h, n + 0-5 by n + 1, and n + 1 by n + 2, these 

equations are immediately recognized in predict or -cor rector terminology to be 
an Euler predictor, an Euler corrector, a Nystrom predictor, and a Milne 
corrector, respectively. Families (l), (2), and (3) are not so accurate as 
the final family and are discarded at the end of each cycle of computation. 

When equations (124) and ( 136 ) are compared, it is clear that 

1 

*i = r 2 = g 


r 3 = 1 

and the coefficients in the difference-differential equations are 


a 2 - “ 

72 = 6 2 = 1 

4 = \ > 

4 = 72 = 0 





* 

7l = 0 , 

72 = 1 


R * _ R * _ 1 o* _ ± 

Bi - S 2 - j , 03 - g 



(137) 
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The coefficients in the operational form are 


m 

1 

2 

bmp 

1 

1 

Rm 

n mp 

-p* 

3 

^m2 

1 

3 

i 

3 

i 

5 

i 

5 

1 

6 

1 

6 


l 

6 


i 

1 



1 

3 

2 

12 



12 


1 




1 

k 

6 




2k 


l 





5 

2? 






(138) 


The error terms are easily calculated. Making use of equation (99 ) > one 
finds for the error in the particular solution 

er ,a = - 5>p 2 + lOA^ - 30 A 3 ] (139) 


and for the error in the complementary solution 


er A = 


24A 5 H 5 

180 


(140) 


Notice that the fourth -order Runge-Kutta method is more accurate as a simple 
integrator (A = 0, er^ = -(pH) 5 /l80) than it is as a differential analyzer 
( \i = 0, er^ = 24(AH) 5 /l80) . On the basis of lowest order error estimates 
(since all methods are referenced to H, they are directly comparable) it is 
not so good as Hamming’s unmodified method (table 11(b)), and a full order 
of magnitude worse than Hamming's modified method (table Il(c)). But, of 
course, when expressed in fundamental families, these are four- and five-step 
methods, respectively. It is comparable in accuracy (46/180 & 0.255) bo a 
three-step, Adams -Bashforth-Moulton, predictor -corrector combination 
(table 11 (b)). 


Stability of Runge-Kutta Methods 

The fourth-order Runge-Kutta approximation to the complementary solution 
of the representative equation can be constructed at once from ( 138 ) and is 

|e - 1 - Ah - | A 2 h 2 - | A 3 h 3 - gq- A 4 h 4 | u n = 0 

There is only one root, the principal one. At first, it might seem that 
stability is not an issue in such a method since there are no spurious roots 
and the principal root approximates e^ h which certainly falls on or inside 
the unit circle for negative A and small enough h. However, if we consider 
the accuracy and stability criterion given by conditions ( 74 b) , the question 
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of stability again arises , even for one-root methods. In fact, in cases for 
which this criterion holds (i.e., when a negative eigenvalue in a set of dif- 
ferential equations has a large absolute value relative to those that are 
driving the system) , the principal root itself may determine the stability 
boundary . 

The variation of the principal root in the vicinity of the unit circle 
is shown in figure 10 for these one-root, one-step methods. In each case the 
points are separated by an increment of 0.1 in h, and |a| = |e iw | = 1. One 
reference value of h is given in each set of points so quantitative esti- 
mates can be made. As indicated, to varies from jt/2 to tz in steps of jt/lO. 

Figure 10(a) shows the principal root behavior for 

Ax = 1 + Ah + | A 2 h 2 

which results from the method formed by combining an Euler predictor with a 
modified Euler corrector. Such a method is self-starting and extremely easy 
to program. It is often used in exploratory numerical research. For real 
negative A the method is seen to be stable for 0 < -Ah = -AH < 2. For 
imaginary A, however, the principal root actually falls outside the unit 
circle for all | Ah | > 0. For imaginary A and h < 0-5 the method is accu- 
rate enough so that more than 200 steps would be required for the instability 
to become serious for most cases. Nevertheless, strictly speaking, its sta- 
bility boundary is zero and, as a completely general method, it should be 
used with caution. (it is unsatisfactory, for example, for studying problems 
containing high-frequency, low-amplitude noise.) 

Figure 10(b) shows the principal root behavior for 

Ax = 1 + Ah + | A 2 h 2 + | A 3 h 3 

which results from third-order Runge-Kutta methods (e.g., Heun's method, see 
ref. page 236) . For real negative A this method is stable for 

0 < -Ah = -1.5AH < 2.5. Now when A is imaginary the principal root falls 
inside the unit circle until h 1.7 • This method has the stability boundary 

| Ah | c = | 1 • 5AH | c = 1.7. 

Finally, figure 10(c) shows the behavior for 

Ai = 1 + Ah + g- A 2 h 2 + ^ A 3 h 3 + — jj- A 4 h 4 


which results from the standard, fourth-order Runge-Kutta method represented 
by equations (136) . This method shows excellent stability for all complex A. 
For real negative A the method is stable for -Ah = -2AH <2.8. The worst 
case occurs when go ^ 0.8rr and limits the stability boundary to 
I Ah | c = | 2 AH | c = 2.6. 


In conclusion, the fourth-order Runge-Kutta method is accurate, easy to 
program, has a higher stability boundary (even when compared on the basis 
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(a) Euler predictor followed by a modified 
Euler corrector. 


( b ) Heun 1 s method . 


2.70 — 



o 

(c) Fourth -order Runge-Kutta method. 
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Figure 10.- Principal root structure of three well-known one-root, one-step methods 



of H) than any of the predictor, one-corrector combinations given in table I, 
and presents no problem in starting or step modification. For many practical 
purposes it is quite satisfactory. 


The Four Iteration, One-Step Incomplete Method in General 

In the simple one-step method the equation for the error in the comple- 
mentary solution, equation (l35)> takes a remarkably simple form. Noting from 
the derivation that (0)° = 1, we see that if er^\ -*• 0(h 5 ), we must have 


1 


-‘02 


(m - 1) 


m = 1, 


5 


( 141 ) 


regardless of the choice of any of the other parameters. These values of 

are the first five coefficients of xJ in the expansion of e x . Hence, 
er^ is always given to the lowest order by 


er A « 


(Ah) 5 24(AH) 5 

5'. lBo 


( 142 ) 


regardless of the choice of the sampling points r 1} r 2 , and r 3 - This proves 
that there exists in equations ( 124) no other one-step, incomplete four- 
iteration method that is more accurate than the standard, fourth-order Runge- 
Kutta method (given by eqs . (136)) in calculating the complementary solution. 
There can be improvements in the accuracy of the particular solution, but, 
since equations (139) and (l40) show the error in the complementary solution 
is nearly the largest, methods that provide these improvements are of limited 
interest . 


Multistep, One-Iteration, Complete Combined Methods 


In this part let us consider the multistep form of the methods described 
above as methods 4a and 4b. These methods are both complete and combined, 
but are easy to analyze and instructive. 


In terms of fundamental families, the difference-differential equations 
for the multi step form of method 4a are 


( 1 ) 

u n+r ' 


k+i 
^ — 1 


L 

j= 2 


a i u n+J + 


a' hu^ 1 ^ ' 
a j nu n+r+i-j 


HR' (1) ’ 

u n+k “ hPiUn+r 


k+i 



Pj u n+J + 


0 


’•hu 


(1) ’ 

n+r+i 



( 143 ) 
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Inserting equation (37) and using operational notation, one can show 


k+i 


k+i 


a =1 


k+x 


E r ' k (E k - Ah I &) - E 


a jB " 


<5=2 


k+i 


-E r " k Ah ^ PjE J E k - f 3 jE J 




U- 


(l) 


n 


Ui 


■n 


= Ahe |ahn e lih(r - k) 


k+i 

E “K hj 


< 3=2 

k+x 


E* 


ij=i 


uhJ 


(i44) 


Compare this with equation (90) • If Oj and Pj are set equal to zero, and 
the bars are deleted from the primed terms, the equations are identical except 
for the term eF” multiplying the left column and the term e ^ k ( r-k ) multi- 
plying the entire right-hand side- The characteristic equations for the two 
methods differ only by the factor E r _k , which means that for given values of 
a and (3, the roots are identical. Further, if we solve only for the final 
family, the factor E r-k has no effect on the particular integral, since it 
appears in both the numerator and denominator of the solution. Finally, the 
jj, h ( r -k ) 

term e simply multiplies the particular integral, so the complete 

solutions to equation (l44) can at once be written from the analysis of 
equation ( 90 ); there results 


2k+x 


E 


ek 


(laj + AhL 2j )E 2k+1 -j 


k+i 


u n = hAe 


qhn 


R 




jxh(k+r+x-j) 


J=2 


i=l 


where 


L xj - Pj > 
R xj = > 


A = 2, k + 1 
a = 1 , k + 1 


(145) 


(146) 


L S j 


1 





j- 1 

^ (Pi a j+i-i 
i=2 


a iPj+3.-i) ’ 


j = 2, 2k + 1 J 


The inverted equations, giving a and 3 in terms of L and R, are the simul- 
taneous, linear equations in (94). 

The only influence of r is contained in the exponent of e in the 
right-hand side of equation (l45) • However, this has a profound effect on 
the results. The equations which must be satisfied for a local polynomial 
fit of order L are now 
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k+i 


Z(k + r + 1 - j) Z_1 Rij + (2k + 1 - j) Z Lij 


1=1 


( 2 k ) 1 = 0 

l “ 0, 1, 


L (147a) 


( 2 k + 1 - orLsj - Rij (k + r + 1 - j) 


= 0 


j- 1 Z = 0, 1, . . L - 1 (147b) 

instead of ( 66 ) and ( 67 ), respectively. The errors are given by evaluating 
the left side of equations ( 147 ) when they are no longer zero and substituting 
the results for the terms inside {} in equations (64) and ( 65 ) . 


To get some idea of the effect of r on the accuracy and stability, one 
can construct tables for the coefficients of L and R for k = 2 and compare 
them with table III. Thus for equations (l47) with k = 2: 


er pi 


l 

Rn 

^12 

Rl 3 

L12 

L13 

( 4 ) 1 

0 

0 

0 

0 

1 

1 

1 

1 

1 

1 

1 

3 

2 

4 

2 

2 ( 2 +r ) 

2 (l+r) 

2 r 

9 

4 

16 

3 

3 ( 2+r ) 2 

3 ( 1+r ) 2 

3 r 2 

27 

8 

64 

4 

4 ( 2 +r) 3 

4 (l+r) 3 

4 r 3 

81 

16 

256 

5 

5(2+r) 4 

5 ( l+r ) 4 

u 

LT\ 

243 

32 1024 

er |J.2 

l 

^22 ^23 

in 

cvi 


E11 

E12 

^13 

0 

1 1 

1 1 


1 

1 

1 

1 

3 2 

1 0 


2+r 

l+r 

r 

2 

9 4 

1 0 

(2+r) 2 

(1+r) 2 

r 2 

3 

00 

kr 

1 0 

(2+r) 3 

(1+r) 3 

r 3 

4 

81 16 

1 0 

(2+r) 4 

(1+r) 4 

r 4 
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Note that this table reduces to portions of table III when r = 2 . A table for 
erx is not constructed, since it is still true that er^ = er^ when A = jjl . 

When k = 2, the characteristic equation for (l44) is 

E — E AI 1 L 22 ) — E^(li3_3 + AhL^s) - Ah[L24E + = 0 (l48) 

and the stability at Ah = 0 is given by the roots to the quadratic 

E 2 - L 12 E - L 13 = 0 (149) 


If we satisfy equations (l47) for L - 4 and various r values between 1.6 
and 2 . 0 , we can construct the curves for L 13 , er^i, and er ^ 2 shown in 
figure 11. The error expressed by er ^ 2 completely dominates the accuracy 

of the method and varies from about 


04- 


( M h ) 5 



\ H 7 h 3 


20 




18 



-0.17V 4 h 5 


20 


at r = 2 to about 
-l.SAq^h 3 at r = 1.6. The varia- 
tion of er^ shown by the solid curve 
is derived using only the first term 
in the truncation error. The actual 
values of er^XlO 5 (for Ah =• 0.1), 
plotted in the circles, indicate that 
the first term in the truncation error 
is accurate for values of Ah = 0.1 
and lower. At around r = I. 635 , Ia 3 
goes to zero and the three spurious 
roots all lie on the origin when 
Ah = 0, giving Adams -Moulton type 
stability. The critical stability 
boundary | Ah | c , see ( 73 ) > is around 
0.3 when L 13 is zero. This boundary 
is determined from plots such as those 
shown in figure 12 . 


Figure 11.- Variation of terv.s controlling stability 
and accuracy of ".ethods given by equations ( 1 U 3 ). 


All things taken into considera- 
tion, the value r = 1.635 appears to 
be a good compromise for this method. 
If we reference the accuracy and stability to H (= 2h) -- according to 
equation (99) --we find 


I am I ~ I (cw) * w ^ 4h= ~ 0 . 0208 V 4 h s 


(150) 


0 < (-AH) C & 0.6 


(15D 


The coefficients in the operational form are, again for r = 1.635 
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Xi 

L mj 

2 

3 

4 

5 

1 

0.98414-6029 

0.01553571 

0 

0 

2 

2.29749376 

-2.44603026 

1.53842844 

-0.37435229 


IXJ 

Rmj 

1 

2 

3 

4 

5 

1 

0.84801929 

0.18240329 

-0 . 01488289 

0 

0 


and using equations (9*0, one finds the corresponding difference-differential 
equations 


(!) 

u n+l . 635 - 


T (l) T T (l) T 

+ a 2 U n+1 + <^ 3 % + O 0 2 hu n+o# 635 + a 3 huh-o .365 


u n +2 

i ( 1 ) 1 

= PihUn+l .635 + p 2 u n+ 

1 + P3U n 

t ( 1 ) ' 

+ P 2 hU n+ o.635 

a R 1 Vi b)’ 

+ P3hU n _o.3 S 5j 


where 

a 2 = -21.44241590 

P 2 = 

0 . 98446029 

(152a) 


a 3 = 22.44241590 

P 3 = 

0-01553571 




a 2 = 20.48107725 

Pi = 

0 . 84801929 

L (152b) 


cl ' 3 = 2.59633865 

P 2 = 

0 . 18240329 





p3 = 

-0.01488289 



A proof that method 4b is unstable for 0(h 5 ) proceeds along the following 
lines. The matrix equation for the operational form of the method, as applied 
to the representative equation (37) > can be written 

r 


1 r 

1 

r 3 



E 3 ** { E 2 — ( cc3+Ahot 2 ) E — Ah.oC/3 } 


-a 2 E 


-E° ~ 2 { Ah| 3 lE 2 + ( p 3 +AhP 2 ) E+AI1P3} E(E-P 2 ) 


u, 


( 1 ) 


■n 


u, 


•n 


=Ahe' aht V ih(r - 2) 


' uh( 3 -j) 


a j e 




LJ 


. 1=1 


pje^ h(3 -j } 
( 153 ) 

One can easily show that the characteristic equation, DE(E) = 0 , reduces to 

E 3 - E 3 (Li2 + AhL 22 ) - E(Lq_ 3 + AhL 2 3) - AhL 2 4 = 0 ( 154 ) 
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where L are simple combinations of the a and p. Now the error in the com- 
plementary solution is determined by substituting e^* 1 for E in equa- 
tion (154) and finding to what order the expanded result does not match the 
expansion of e^ 1 itself. This simply amounts to entering table IV for 
k = 3 and finding the first row that does not sum to zero. If we are to make 
the first five rows all sum to zero, so the order of the error is h 5 , all 
five values of L are completely determined; they must satisfy the equations 


Li 2 - -8 

L 13 = 9 

L 22 = 17/3 
L 2 3 = lV3 
L 2 4 = -1/3 


(155) 


regardless of the values chosen for r and R. Using equations ( 155 ) and 
setting h = 0, we find equation (15^-) reduces to 


E 2 + 8E - 9 = (E - 1)(E + 9) = 0 


which has a violent instability given by the root E = - 9 . 


Two -St ep, Two -It erat ion, Incomplet e , 

Combined Methods 

The final section of this part is devoted to a discussion of methods 2a 
and 2b. Method 2a is really a special form of equations (12^) when the latter 
are simplified to two iterations, or one value of r- Applying the analysis 
of those equations, one can show that all methods of type 2a having an error 
of 0(h 5 ) have operational coefficients given by (the * has been omitted from 
R*i in the following) 


v 

L mj 


R mj 


m\ 

2 

3 

1 

2 

3 

*1 

8(r - 2) 

17 - 10r 

2 

2(ij-r 2 -9 r + 4) 

-2(2r 2 - 4r + 1) 

X 

1 - 2r 

1 - 2r 

r ( 1 - r ) ( 1 - 2r ) 

(l - r) (l - 2r) 

u 

CJ 

l 

1 — 1 
U 

O 

-4(r - 2) 

2(5 - 4r) 

0 

-2r 

2(1 - r) 

d. 

1 - 2r 

1 - 2r 

1 - 2r 

1 - 2r 

O 

-2r 

2(1 - r) 




O 

1 - 2r 

1 - 2r 





and difference -differential equations given by 
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u n+r _ a 2 u n+i + a 3% + h(a 2 u n +i. + 0-3 Un.) 


f 1 ) ' 

u n+2 = P2 u n+i + 33U-n + h ( Piu n + r + P 2 u n+i + P3 u n 


(156a) 


where 


a 2 = r 2 (3 - 2r) P 2 = 8(r - 2)/(l - 2r) 

a 3 = (l + 2r) (l - r) 2 p 3 = (17 - 10r)/(l - 2r) 

Pi = 2/[r(l - r) (l - 2r) ] 

P2 = 2[4r 2 - 9r + 4]/[(l - r)(l - 2r) ] 
P 3 = -2[2r 2 - Ur + l]/[r(l - 2r) ] 


a-2 = -r 2 (l - r) 
0,3 = r(l - r) 2 


(156b) 


The leading error terms are 

[ ( 5r 2 - llr + 4 )\i + 5r(l - r)A]p 4 H5 


erp. = 


er A = 


l80[2(2r - 3)] 
(2 - 3r)A 5 H 5 


180 (2r - 3) 

and the characteristic equation at h = 0 is given by 


(157) 


( e - d(e + : g° r ) = 0 


( 158 ) 


so the spurious root starts at 

A 2 


17 - lOr 
1 - 2r 


which vanishes at r = 1.7, giving Adams -Moulton type stability there. A 
table of the error terms for various values of r is included and stability 


r 

erp.1 

pSRS 

erp.2 

Ap 4 H 5 

er A 

A 5 H 5 

2.00 

0.0056 

-O.O278 

-0.0222 

1.90 

.0040 

-.0297 

-.0257 

1.85 

.0030 " 

-.0312 

-.0282 

1.80 

.0019 ~ 

" -^0333" 

-.0314 

1-75 

1.70 

.0003 

-.0364 

-.0361 

- . 0017 

-.0413 

-.0430 
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Method 2b can be written 



and reduced to the operational form 

{E 3 - E 2 (L 12 + AhL 22 + A 2 h 2 L 32 ) - E(AhL 23 + A 2 h 2 L 33 ) - A 2 h 2 L 34 }u n 


(159) 


— Ahe* 1 [ (Rh + AhR 21 ) e 2 !- 1 ' 1 + (Ri 2 "*■ AhR 22 )e^ 

+ AhR 23 + R 11 e |jJ:1 ^ r+1 ^ + R^e 111 ^] (l60) 


where 

Rn = ft 2 

^12 “ P3 

Rix = Fl 
R 12 = P 2 

R 2 i - L 32 

R22 = L33 

R 2 3 ~ L 34 


L 12 - ^2 

-^22 ~ 3 2 ^2 ^ 2^1 

L32 ~ a 2 Pl 

1^23 “ 33 “ 32 a 2 + 32 a 2 
L 33 = -P 2 a 2 + Fi a 3 + P 2 a 2 
Ii34 = P 2 a 3 - P 3 a 2 


(l6l) 


The accuracy conditions formed by collecting the terms independent of Ah in 
equation (l6o) are 


er pi = 


(ph) 

V. 


RijU3- J) Z " 1 +R 1 jUr+2- j) 


Z-1 


+ 2 l L 12 - 3 l 


l — 0, 1, 




k+i 


(l 62 a) 


The term corresponding to ^(j-l)Lij in equations (64) and ( 65 ) is unity. 

The second set of conditions formed when the coefficients to Ah are 
collected in equation ( 166 ) are 
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er M o = 


Ahpih) 


M - 2 Z! 


2 

R 2 jl(3 - j) Z_1 - ^ [Rij(3 


1=1 


d) ^ + Ri j( r + 2 - j) ] 


1=1 


3 

y - j: 


1=2 ■ 


Z = 0, 1, . . . 


(162b) 


The coefficient to (Ah) in equation ( 60 ) is identically zero. 


If erp. is made to be of order h 5 , then equations (162) provide nine 
equations for the ten unknown constants (four a,, five | 3 , and r) in equa- 
tions ( 159 ). Hence, one can calculate the errors er^ and er ^ 2 as func- 
tions of the parameter r. The result is shown in figure 14 , where, Twe see 
that erp.2 completely dominates the total error except in the region r~1.8. 



. 6 - 

4- 

|XH| C 

2 - 



1.5 


1.7 


1.8 


1 9 
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Figure 1*4-.- Variation of terms controlling stability 
and accuracy of methods given by equations (159) • 



Notice, from equation (l6o) , that the method has Adams - Moulton type stability 
for all r, since at h = 0 the characteristic equation reduces to E(E - 1) = 0. 
If r is set equal to 1.8, we find the following values for L and R 


v 

L mj 

R mj 

R mj 

"\ 

2 

3 

4 

1 

2 

3 

1 

2 

1 

720 
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0 

705 

15 

0 

375 

-375 

2 

-1128 

1848 

0 

1340 

932 

-64 

0 

0 

3 

1340 

932 

-64 







Divide by 720 


er. . = 


(laH ) 5 , er A = ^ ( M ) 5 


L M- 720 720 

and the corresponding set of difference-differential equations 


Un+i.s _ u n+i + (^ 268 u n+ i + 22u n - 230u^+o.8 


u n+2 - u n+i + Jjg (Wu n +i + + 25 un+i.s - 25 u n+0 .8 


( 163 ) 


The error terms for these equations are by far the smallest errors 
(l/720 0.0014) for any of the methods considered above. However, there is 

the usual sacrifice in stability. From a study of results such as those 
given in figure 15, one can calculate the curve for the induced stability 
boundary shown in figure 14. We see that equations ( 163 ) are limited by the 
stability boundary |AH| C - 0.3* 


THE OPERATIONAL FORM 


Definition and Discussion 

A variety of systems of difference -differential equations have been 
analyzed as they applied to the solution of ordinary differential equations. 

In every case the method being studied was associated with an equation of the 
form 

P(/mjPh) m_1 E k+1 ”5) = Ahe^ hn Q(R m j(Ah) m-1 e |J ' h ^ ri+1 "^^ (l64) 

where P and Q symbolize polynomials with terms such as those within the 
arguments. In every case this equation was the sole basis for determining the 
accuracy and stability of the method. In general, the maximum value of m 
is determined by the number of iterations. If a method uses M iterations. 
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m max = M + 1. The maximum value of j is k + 1, where the interpretation of 
k varies according to the type of method. For incomplete, uncombined meth- 
ods, such as Hamming* s, Adams -Moulton, etc., k is the largest step number 
used in the predictor or the corrector. In complete, uncombined methods, k 
is twice the maximum step number. In combined methods, steps are not neces- 
sarily equispaced and k loses its connection with step number. The term rq 
can be replaced with k in the uncombined methods. In the combined ones, 
however, rq determines the location at which calculations of the function 
and/or its derivative are carried out and these locations are not necessarily 
spaced in integer multiples of h. 

We refer to equation (164) as the operational form of a numerical method 
and to the coefficients and L m j as the coefficients in the operational 

form. One can show: 


1. All linear, dif fere nee -differential equations with constant coeffi- 
cients have an operational form. 

2. Any two such methods with the same operational form have, except for 
round-off considerations, the same accuracy and stability and give, 
therefore, except for round-off considerations, identical numerical 
results when applied to equations (ll) . 

3- The coefficients 7, • • • in the actual difference- 

differential equations affect the accuracy and stability of the 
method only as they affect the coefficients in the operational form 
and the correspondence is not unique. 


Let us consider in more detail the above statement number 2. Strictly 
speaking, two methods that give identical numerical results when applied to 
linear equations cannot be classified as equivalent, since one can require 
more iterations than the other. As an example, the incomplete, four-step 
method. 


u n +4 - 3-9 u n+3 “ 5-7 u n+2 + 3-7 u n+i - 0-9 u n 


u n+4 - Un+3 + 


J7T ( 5 u n+4 


8 U ; 


•n+3 


“ Un+e") 


(165) 


when applied to equations (ll) gives results 15 that are identical with those 
of the complete two-step method presented in ( 96 ) . Notice that the number of 
multiplications and storage requirements of the two methods are the same. 
However, the use of equations ( 165 ) requires twice the number of iterations 
and for this reason, as a numerical method, is neither equivalent nor 
practical. 


15 To show this numerically, care must be taken with the initial 
conditions since neither method is self-starting. 
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II 


Accuracy 

The accuracy which a given operational form provides in the solution of 
a set of linear differential equations can be measured by the magnitude of the 
error terms er^,, defined by equations (63) and (57); and by the error term 
er^ where er^ = (er^)^^. We are concerned only with polynomial approxi- 
mation, so these terms are expanded in powers of h and the lowest power of 
h with a nonvanishing coefficient gives the order of the local polynomial 
fit. 


Clearly, the parameters \x and A depend on the differential equation, 
so the error of any method must be expressed as a function of jjl, A, and h. 
First, we notice that the error in the particular solution can always be 
expressed as a polynomial in Ah. Collecting terms in this way, we next see 
that the coefficients of this polynomial are functions of ph and (h-* 1 ) . 
Finally, these coefficients are each expanded in powers of ph and made to 
vanish to the desired order. This leads to sets of conditions on the coef- 
ficients in the operational form, examples of which are given in 
equations (66), (67), and (131) through (13^). 


Stability 

In this part we show that the Dahlquist stability theorem, derived for 
an implicit multistep equation, also holds for any multistep, predictor - 
corrector method contained in equation (l64), provided r± = k. That is, 
provided equispaced , predictor - corrector methods are not combined with Runge - 
Kutta techniques . 

The argument starts by inspecting equations (3) and (A). As we have 
seen, the degree of the polynomial embedded in equation (3) depends upon how 
many of the terms erp(O), er p (l), . . . , erp(L) in equation (A) can be made 
identically zero. If, in fact, the P values are chosen so all terms 
through erp(L) are zero, L + 1 equations or constraints on (3 must be satis- 
fied and the order of the embedded polynomial is L. On this basis, the 
maximum value of L is twice the step number. The stability, on the other 
hand, depends on the roots to the characteristic equation derived from 
equation (3 ) , namely, 

k+i 

E k - ^ (Pj + Ah|3j)E J = 0 (166) 

j=i 

The Dahlquist theorem is based on the hypothesis that, for a method to be 
stable, it is necessary that the roots to the characteristic equation lie on 
or inside the unit circle in a complex plane when h = 0. In our notation 
and analysis, these statements mean that L is to be made as large as 
possible in 
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I 


k+i 

Y, z ( k + 1 - O) 2 _1 Pj + (k + 1 - 3 ) Z Pj = k Z , * » 0 , 1 , . . ., L (l6?a) 

3 =1 

without having the absolute value of any of the roots to 

k+i 

E k - y | 3 jE k+1 "' 3 = 0 (167b) 

3=2 

exceed unity. The Dahlquist theorem states that there are no combinations of 
(B for which this is possible if L > k + 2 (k even) or L > k + 1 (k odd) . 

The extension of these results to uncombined methods constructed from an 
operational form with Tj_ = k is quite simple. Consider, for example, the 
two-corrector method expressed by equations (ill) . The error, er^, for such 
a method is given by equations (117). For h = 0 we find a necessary condi- 
tion for obtaining a method embedding a polynomial of order L is 

k+i 

y Z(k + 1 - j) Z_X Ri j + (k + 1 - j) Z Lo.j = k z , l = 0 , 1 , . . L (l68a) 

and a necessary condition for stability is that the absolute value of the 
roots to 

k+i 

E k - y LijE k+1-J = 0 (l68b) 

3=2 

does not exceed unity. 

Examining the discussion of equations (167), we see that the Dahlquist 
stability criterion is quite independent of the role of (3 in the difference- 
differential equations. Thus, although L and R of equations (l68) can enter 
the difference -differential equations entirely differently, the Dahlquist 
theorem immediately tells us that only the first k + 2 (for even k, or k + 1 
for odd k) of equations (168a) can be satisfied if the absolute value of the 
roots to equation (l68b) is to be no greater than one. Hence, the Dahlquist 
theorem still applies to combined predictor -corrector methods when any number 
of correctors, all of which may be different, are used in an uncombined, 
multistep, predictor -corrector method. 

If combined predictor -corrector formulas are used (complete or incom- 
plete), r-f is no longer equal to k. For example, equations ( 12 ^) have, in 
place of (l68a), the accuracy conditions 
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k+i 


Z(k + 1- j) Z-1 L 1 j+ (k + l-j)Vj 




Z -i * 


Rii=k 


J=2 


X=1 


l = 0 , 1 , 


(169) 


although the stability equation remains identical to (l68b) . Using the same 
argument as above, we see that the Dahlquist theorem is no longer applicable 
to these cases. 

The fact that the Dahlquist stability criterion must be modified if 
unequal steps are taken in advancing the integration of differential equations 
has been recorded by several authors, for example, references 13 through lb. 

But the variety of meanings given to the words "step number” in these and 
other references complicates a comparison of the stability capabilities of 
the various methods in a sense similar to that studied by Dahlquist. This 
problem, already discussed, is discussed here in light of its connection with 
the Dahlquist theorem. 

Consider first the "conventional," four -step, method composed of an Adams - 
Bashforth predictor (line 5> table l(a)), followed by an Adams -Moulton correc- 
tor (line 4, table l(b)) as symbolized in sketch (g) . (The symbols are 
defined in the previous section.) 




Sketch (g) 


The data used to calculate the value of the predicted and corrected function 
are encircled; the remaining data are ignored. The step size is shown as h a , 
the choice which coincides with the definition (123) • It is equal to H, the 
distance advanced by two iterations. Six bits of data are weighted in the 
corrector and a stable method results with a local polynomial fit of order 
five. The error (see table Il(j)) is around -0 . 

Consider next the combined, two-step method presented by Butcher in 
reference 15* In our notation it reads 
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u n+i .5 = u n + g h b (9un+x + 3u„) 


u n +2 - ^ ( 2 8u n+1 - 23 u n ) + h b ^32u^ + i #5 - 60^+! - 26u n 


,(i>’ 


> (170) 


u n +2 = (32u n+i 


> 1 , 4 , C 1 )' , c ( 2 )' , 0 1 

u n ) + — h b ( 64 u n+1 . 5 + 15 u m . 2 + 12u n+ i 


un 


where also coincides with the step size defined by ( 123 ). The process 

is symbolized in sketch (h) • The error of the method is around (Ahb) 6 /l 24 





or ~ 0 . 06 (AH) s , since hb - 3 H/ 2 . The increase in accuracy is compensated 
for, as usual, by a decrease in the stability boundary. 


Although the method illustrated by sketch (h) is by definition a two- 
step method, it appears very similar to the four -step method shown in 
sketch (g). The real connection between the two is discovered by reducing 
each to its operational form. For simplicity, re -reference the indexing in 
equations (170) to a step size equal to h a so that, for example, the first 
equation reads 


0 

u n +3 - u n 


^a 

16 


(9%i+ 2 + 3u a ) 


When this is done one can show that both methods are governed by equa- 
tions (l68) in which k - 4 - Apply Dahlquist T s theorem to the latter, and we 
find that a local polynomial fit of order 5 can certainly 16 be found which 

ls Actually the theorem states that a local polynomial fit of order 6 
would be stable, but such methods would have a spurious root on the unit 
circle at h = 0, and are usually unstable for h > 0. 
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is stable for a variety of methods containing, over the length 2hfr, the kind 
of data shown in sketch (h) . Butcher showed that the particular choice of 
data encircled in sketch (h) are, in fact , stable for a polynomial of order 
5 when weighted as in equations (17O). If the number of steps, sized hfc, is 
increased to 3, a method similar to that given by equations (17O) — in that 
the values of the function and its derivative at the midpoints of steps 
behind the last are ignored -- is stable according to Butcher for a polynomial 
of order 7. However, Butcher also showed that when kb, the number of h^ 
steps, is increased beyond three, this process is unstable for polynomials of 
order 2kb + 1. On the other hand, if we use the operational form, the 
Dahlquist theorem tells us that some choice of such equispaced data in the 
same interval can be used to derive a method that is_ stable for 2k^ + 1 or 
even 2kb + 2 in the sense described in footnote l6. 

Of course, the above example, when used in this perspective, is quite 
unfair as a true measure of the value of combined methods. This is due to 
the symmetrical choice of the sampling points. But it is quite useful in 
demonstrating that the words "step-number , " as they are used in contemporary 
literature, are to be treated with great caution in comparing methods. 

Perhaps the simplest way to present the basic issue involved is to con- 
sider sketch (i) . The implicit method defined advances the solution a step 

h in one cycle of computation and, on this basis, 
it is a two-step method. The generalization of 
Dahlquist ! s stability theorem will provide an 
answer to the following question: 

Is there any value of r for which the 
data encircled in sketch (i) can be used 
to construct a stable method embedding 
a local polynomial of order 8? 

Nine bits of data are now available so that a 
method having a polynomial fit of order 8 can 
easily be constructed. The Dahlquist theorem, when applied to the opera- 
tional form, immediately tells us that the answer to the above question is 
negative when r = 3/2. If r = 4/3, however, the situation is not so simple. 
Then the data are spaced so as to be identical with operational forms of 
certain equispaced methods with six steps. The Dahlquist theorem tells us 
that there are some operational forms of such cases that are stable, but 
whether or not they can be obtained omitting data at two of the intervals is 
not known. Admittedly, questions such as this are approaching the academic, 
but their answer is a fundamental aspect to one area of numerical analysis. 



Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 94035, Jan. 24, 1967 
129-04-03-02-00-21 
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TABLE I.- COEFFICIENTS IN DIFFERENCE -DIFFERENTIAL EQUATIONS FOR CERTAIN PREDICTOR -CORRECTOR FORMULAS 


(a) Predictor formulas 



a*. 

oa 

a 2 

1 

a 2 

a 3 

a 3 

a 4 

a; 

a 5 

cx 4 

Os 

a p 

1 

Euler 



1 

1 








l/ 2 h 2 u lx 

2 

Nystrom 




2 

1 







l/ 3 h 3 u llx 

3 

A-B* two-step 



1 

3/2 


-1/2 






5 /l 2 h 3 u l:L1 

4 

A-B* three -step 



1 

23/12 


-16/12 


5/12 




3 / 8 h 4 u iv 

5 

A-B* four -step 



1 

55/24 


- 59/24 


37/24 


- 9/24 


25 l/ 720 h 5 u V 

6 

Milne -Hamming 
(No mod.) 




8/3 


- 4/3 


8/3 

1 



i 4 / 45 h 5 u v 

7 

Hamming (mod.) 



1 

8/3 


-4 


4 

1 

- 8/3 

-1 

-ll/72h s u y 

8 

Crane -K 



1.547652 

2.002247 

-1.867503 

-2.031690 

2.017204 

1.818609 

- 0.697 353 

-0.714320 


0 . 40 l 6 h 5 u y 

9 

( 0 . 01 , - 0 . 01 ) 



-59/17 

127/34 

76/17 

59/34 






3 / 68 h 3 u ll:L 

10 

I Stetter 



-4 

4 

5 

2 






l/ 6 h 4 u iv 


*Adams-Bashforth method 


(b) Corrector formulas 



Pi 

Pi 

02 $2 

33 

Pi 

P 4 

Pi 

Ps . Pi 

Pe 

Pp 

1 

Modified Euler 


1/2 

1 

1/2 








-l/l2h 3 u 1:Li 

2 

A-M* two-step 


5/12 

1 

8/12 


-1/12 






-l/ 24 h 4 u^ 

3 

A-M* three-step 


9/24 

1 

19/24 


-5/24 


1/24 




-19/720h 5 u v 

4 

A-M* four -step 


251/720 

1 

646/720 


-264/720 


106/720 


-19/720 


-3/l60h s u v:L 

5 

Milne 


1/3 


4/3 

1 

1/3 






-l/90h 5 u v 

6 

Hamming 
(No mod.) 


3/8 

9/8 

6/8 


-3/8 

-1/8 





-l/ 40 h s u v 

7 

Hamming (mod.) 


42/121 

126/121 

108/121 

0 

-54/121 

-14/121 

24/121 

9/121 



-2l/l210h e u vi 

8 

(0.01, -0.01) 


34/93 

12/31 

100/93 

19/31 








9 



3/ll 

-27/11 



27/11 

1 

3/11 






* A dams -Moulton method 





TABLE II. - COEFFICIENTS IN THE OPERATIONAL FORM OF A NUMBER OF METHODS 

(a) Predictor, row 6 of table l(a) (Milne-Hamming (no mod.)); 
corrector, row 5 of table l(b) (Milne) 


\j 

L mj 



R mj 


m \ 

2 

3 

4 

5 

6 

l 

2 

3 

4 

5 

6 

i 

0 

9 

0 

0 


3 

12 

3 

0 

0 


2 

12 

3 

0 

3 


0 

8 

-4 

8 

0 


3 

8 

-4 

IT 

0 









| Divide by 9 | 

er^ = (0.0056 - 0.052AH) (|jlH) 5 
er A = 0.0056(AH ) 5 


(b) Predictor, row 6 of table l(a) (Milne-Hamming (no mod.)); 
corrector, row 6 of table l(b) (Hamming (no mod.)) 


\J 

Lmj 

R mj 


2 

3 

k 

5 

6 

1 

2 

3 

4 

5 

6 

1 

9 

0 

-1 

0 


3 

6 

-3 

0 

0 


2 

6 

-3 

0 

3 


0 

8 

-4 

8' 

0 


3 

8 

-k 

8 

0 









Divide by 8 


erp. = (O.O 33 - 0.155AH) (tiH ) 5 
er A = 0.033(AH ) 5 


(c) Predictor, row 7 of table l(a) (Hamming (mod.)); 
corrector, row 7 of table l(b) (Hamming (mod.)) 


\ j 

Lmj 

R mj 

m \ 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

1 

126 

0 

-14 

9 

0 

42 

108 

-54 

24 

0 

0 

2 

150 

-54 

24 

42 

-42 

0 

112 

-168 

168 

-112 

0 

3 

112 

-l68 

168 

-112 

0 








Divide by 121 


er^ = (0.0175 - 0.109 AH)(m-H ) 6 
er A = 0.0175(AH ) 6 
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TABLE II. - COEFFICIENTS IN THE OPERATIONAL FORM OF A NUMBER OF METHODS - 

Continued 

(d) Predictor, row 3 of table l(a) (A-B two step); 
corrector, row 2 of table l(b) (A-M two step) 


\ 3 


L 

mj 

| 

m \ 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

1 

24 

0 




10 

1 6 

-2 




2 

2 6 

-2 




0 

15 

-5 




3 

15 

-5 











Divide by 2k 


er^ = (0.042p - 0.174A)p 3 H 4 
er A - -0.132(AH ) 4 

(e) Predictor, row 2 of table l(a) (Nystrom) ; 
corrector, row 2 of table l(b) (A-M two step) 


\ i 

L tnj 

R mj 

m \ 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

1 

12 

0 




5 

8 

-1 




2 

8 

4 




0 

10 

0 




3 

10 

0 











Divide by 12 


er^ = (0.043-i - 0.138A)p 3 H 4 
er A = -0.097 (AH ) 4 

(f) Predictor, row 9 of table l(a) ; 
corrector, row 8 of table l(b) 


\ 3 


R mj 

m \ 

2 

3 

4 

5 

6 

1 

2 

3 

4 ] 

5 

6 

1 

1224 

1938 




1156 

3400 

544 




2 

-612 

5712 




0 

4318 

2006 




3 

4318 

2006 











Divide by 31^*2 


er^ = 0.0l(p - A)p 3 H 4 
er A = -0.027 (AH ) 5 
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TABLE II • - COEFFICIENTS IN THE OPERATIONAL FORM OF A NUMBER OF METHODS 

Continued 

(g) Predictor, row 4 of table l(a) (A-B three step); 
corrector, row 2 of table l(b) (A-M two step) 


V 


L 

mi ___ 



R 


A 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

l 

144 

0 

0 



60 

96 

-12 

0 



2 

156 

-12 

0 



0 

115 

-80 

25 



3 

115 

-80 

25 










Divide by 144 


er^ = (0.042 - 0. 156AH) (pH) 4 
er A = 0.042(AH) 4 


(h) Predictor, row 4 of table l(a) (A-B three step); 
corrector, row 3 of table l(b) (A-M three step) 


\j 

L mJ 

R mj 

m \ 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

1 

288 

0 

0 



108 

228 

-60 

12 



2 

336 

1 

0 

12 



0 

207 

-144 

45 



3 

207 

-144 

45 










Divide by 288 


er = (0.026p - 0.l4lA)p 4 H 5 
er A = -0. Il4( AH) 5 

(i) Predictor, row 5 of table l(a) (A-B four step); 
corrector, row 3 of table l(b) (A-M three step) 



Lmj 

R mj 

m \ 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

1 

192 

0 

0 

0 


.72 

152 

-40 

8 

0 


2 

224 

-40 

8 

0 


0 

165 

-177 

ill 

-27 


3 

165 

-177 

111 

-27 









Divide by 192 


er^ = (0.026 - 0.131AH) (pH) 5 
ev\ = 0.026(AH) 5 
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TABLE II. - COEFFICIENTS IN THE OPERATIONAL FORM OF A NUMBER OF METHODS - 

Concluded 

(j) Predictor, row 5 of table l(a) (A-B four step); 
corrector, row 4 of table l(b) (A-M four step) 



Lmj 

Rmj 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

1 

17280 

0 

0 

0 


6024 

15504 

-6336 

2544 

-456 


2 

21528 

-6336 

2544 

-456 


0 

13805 

-14809 

9287 

-2259 


3 

13805 

-14809 

9287 

-2259 









Divide by 17280 


er^ = (0.019 m, - 0.122A)m 5 H s 
er A = -0.103(AH ) 6 


(k) Predictor, row 10 of table l(a) (Stetter); 
corrector, row 5 of table l(b) (Milne) 


\ 3 

Lraj 

Rmj 

m \ 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

1 

0 

3 




1 

4 

1 




2 

0 

6 




0 

4 

2 




3 

4 

2 











Divide by 3 


er^ = (0.0056 m - 0.028A)m 4 H 5 
er-^ = -0.022(AH) 5 
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TABLE III. - COEFFICIENTS OF L AED 



k = 1 






2 




R11 


3 



R11 

Rl 2 


4 


R11 

r 12 

Rl3 

l 

5 

Rn 

R12 

Rl3 

Hl4 

0 


0 

0 

0 


1 


l 




2 


10 



4 

a 


75 

48 

27 

12 

i 


500 

256 

1 108 

32 

i 


3125 

1280 

l 405 

: 80 

r«i 


18750 

1 6l44 

1458 

192 































































TABLE III. - COEFFICIENTS OF L AND R FOR USE IN THE CALCULATION OF er^ FOR ONE- THROUGH FIVE-STEP 

METHOD - Concluded 

(b) Equation (67) 



5 : 

-R11 : 

-Rl 2 ;; “Ri 3 

-R3.4 -Rl 5 


.3 ; ' ^22 

R 2 3 1 

^22 

1*23 


Rl3 “Rl4 


Rl 3 • -Rl 4 -Rl 5 


' R22 i ^23 R24 


R pp R23 R24 R25 


R 2 3 R24 R25 R 


L22 l 2 3 l 2 4 


L pp L 2 3 • 1*24 L 2 5 


L 22 1*23 L 24 L 2 5 L 2 s 


2 1 


625 256 8l 16 


3125 1024 243 32 


15625 | 4096 I 729 I 64 


6 4 2 


256 I 108 32 4 


1280 405 80 5 


6i44 1458 192 6 


1 

0 

1 

0 

1 

0 

1 

0 

1 

0 
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TABLE IV. - COEFFICIENTS OF L FOB USE IN THE CALCULATION OF er A ONE- THROUGH FIVE -STEP METHOD 

(see eq. (72)) 



k= 1 



HI 







■ 



L12 

L 2 2 



2 










Ll2 

l 22 

^32 

Ll3 

^23 

L33 

3 







L12 

" 1 ■ 

L 2 2 


Li3 

L23 

L33 

Li4 

^24 

L34 

1 + 




Lx 2 

1 

L 2 2 

^32 

Ll3 

L23 

L33 

Li4 

L^4 

L34 

LlS 

L25 

L35 

k z 

l 

5 

Ll2 

L 2 2 

L32 

Lis 

^23 

L33 

Lj_4 

Lg4 

L34 

Ll5 

L25 

L35 

Lie 

L 2 6 

L36 

k»l 

2 

3 

4 

5 

0 


1 

0 

0 

1 

0 

0 

1 

0 

0 

1 

0 

0 

1 

0 

O' 

1 

1 

1 

l 

1 

1 


4 

1 

0 

3 

1 

0 

2 

1 

0 

1 

1 

0 

0 

1 

0 

1 

2 

3 

4 

5 

2 


16 

8 

2 

9 

6 

2 

4 

4 

2 

1 

2 

2 

0 

0 

2 

1 

4 

9 

16 

25 | 

3 


64 

48 

24 

27 

27 

18 

8 

12 

12 

1 

3 

6 

0 

0 

0 

1 

8 

: 27 

64 

125 : 

4 


256 

256 

192 

81 

108 

108 

16 

32 

48 

1 

4 

12 

0 

0 

0 

1 

16 

81 

256 

: 625 1 

5 


1024 

1280 

1280 

243 

405 

540 

32 

80 

160 

1 

5 

20 

0 

0 

: 0 

1 

32 

243 

1024 

3125 

6 


O 

vo 

OV 

-4- 

-=J- 

1 — 1 

VO 

->] 

OV 

00 

O 

729 

! 

1458 

2430 

64 

192 

480 

1 

6 

30 

0 

0 

0 

1 

1 

64 

729 

4096 ' 15625 






"The aeronautical and space activities of the United States shall he 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof .” 

— National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and technical information considered 
important, complete, and a lasting contribution to existing knowledge. 

TECHNICAL NOTES: Information less broad in scope but nevertheless of 
importance as a contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: Information receiving limited distribu- 
tion because of preliminary data, security classification, or other reasons. 

CONTRACTOR REPORTS: Scientific and technical information generated 

under a NASA contract or grant and considered an important contribution to 
existing knowledge. 

TECHNICAL TRANSLATIONS: Information published in a foreign 
language considered to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information derived from or of value to NASA 
activities. Publications include conference proceedings, monographs, data 
compilations, handbooks, sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION PUBLICATIONS: Information on tech- 
nology used by NASA that may be of particular interest in commercial and other 
non-aerospace applications. Publications include Tech Briefs, Technology 
Utilization Reports and Notes, and Technology Surveys. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION DIVISION 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 



